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INTRODUCTION 

The GSFC Faculty Fellowship Program is a cooperative effort between 
Goddard, the American Society of Engineering Education (ASEE), the University 
of Maryland, and the Catholic University of America. 

Under this program university faculty members come to Goddard for a 
period of ten weeks during which they engage actively in research and at the 
same time participate in seminars related to their research. 

The objectives of this program are many; some of them are listed below: 

1 . Stimulation of schools to become interested in the research problems 
confronting Goddard. 

2. Creation of interest on the part of the university faculty to continue their 
research after completing the formal program. 

3. Stimulation of our people professionally through associations with the 
university faculty and through participation in the program’s seminars. 

4. Establishment of closer ties with the Universities. 

This report marks the completion of the third year of this program at 
Goddard, and represents the progress in one of the objectives of the Summer 
Faculty Program. 
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A SURVEY ON SOME MATHEMATICAL MODELS OF THE 
VERY LOW FREQUENCY WAVE PROPAGATION 

by 

David H. S. Cheng r 


INTRODUCTION 


Historical Background 


The development of the theory of the very low frequency wave propagation 
may be divided roughly into two periods. The first period was ushered in with 
the introduction of long distance communications by means of the wireless at the 
turn of the century. The underlying principle in almost all the earlier analyses 
was the diffraction concept advanced by Rayleigh in 1904 [1]. 

In such a diffraction model, the existence of the ionosphere was ignored, and 
the earth was considered as a homogeneous conducting sphere. Among the earlier 
investigators were Poincare [2], MacDonald [3], Nicholson [4], Sommerfeld 
[5] , Love [6], March [7] , Rybczynski [ 8] , and others. Most of the theoretical 
investigations made use of the theory of the zonal harmonics, either by a direct 
summation or by forming an integral from the series , to explain the facts con- 
cerning the propagation of radio waves around the earth’s curvature. However, 
the series then formulated was difficult to sum because of its slow convergence 
for the radio problem. 

Poincare [ 2 ] and Nicholson [4] attempted to replace the series by an inte- 
gral and then obtained an approximate value for the integral by means of the 
calculus of residues. Their scheme, though substantially sound, was exceedingly 
elaborate, and their approximations became invalid in the vicinity of the antipodes 
of the transmitter. 

The method used by MacDonald [ 3 ] was to approximate the original series 
by a modified series which in turn was replaced by an integral. This procedure 
was questioned from a purely mathematical point of view, though worked in most 
practical circumstances. 

Love [ 6 ] undertook the laborious task of evaluating the series for the radio 
problem by numerically adding the terms of the series taken in groups. His 
memoir contained a complete bibliography of the problem up to 1915. 


^University of Missouri, Columbia, Missouri 
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The results of the aforementioned authors did not explain the propagation 
phenomenon round the earth’s curvature; nevertheless, a correct conclusion 
was reached that the diffracted field decayed exponentially from the transmitter. 

Sommerfeld [ 5 ] investigated the problem using a planar model and con- 
cluded that the field decayed approximately as 1/Vcf or at most 1/ d where d is 
the distance of propagation from the transmitter, which checked fairly well with 
the observed facts. Upon his suggestion, March [7] reformulated the problem 
using a spherical model. He cast his solution in an integral in order to obtain 
an answer more consistent with observed facts. Indeed, he found that the field 
decayed approximately as 1/^6 sin 6 ~ l/d where 9 is the angular distance be- 
tween the transmitter and the receiver. March concluded that the waves were 
propagated without exponential attenuation contrary to the results obtained by 
other investigators. His method was severely questioned. However, von Rybczynski 
[ 8] extended March’s method and pointed out that the exponential factor was 
neglected by March. 

In reference to these investigations, Smith [9] remarked that "many years 
of work of the most famous mathematicians were inadequate to derive this task, 
which at the first glance seems very simple, from that state with respect to 
which Nicholson said that out of all mathematical problems it is the only one 
that has caused such a great divergence of opinions." It was against this back- 
ground that Watson [ 10] , in 1918, succeeded in his first of two classical papers 
on VLF wave propagation in making the final conclusion that we could not explain 
the fact of the long distance propagation of long waves around the terrestrial 
sphere using only the diffraction model. 

The second period was initiated in 1919 by the work of Watson [10]. In the 
second of his two papers on the radio problem, Watson postulated a concentric 
plasma of high conductivity to account for the ionosphere which is mainly re- 
sponsible for the VLF propagation over great distances. Watson succeeded in 
reformulating the integral used by others by means of now the well known Watson 
transformation. Instead of working with the slowly convergent series he was 
able to replace it by a rapidly converging series. He introduced the earth iono- 
sphere waveguide concept. The radio waves excited by a vertical Hertzian dipole 
propagated in the atmosphere between the earth and the ionosphere in a way 
similar to the microwave in modern wave guides . 

Among the earlier investigators in the second period in the development of 
the theory of the propagation of LF and VLF waves, aside from Watson, are 
Rydbeck[ll] , Van der Pol and Bremmer [12], Eckersley [13 ] , Bremmer [14], 
Kendrick [15], Weyrick [16 ], Brekhovskikh [17], and Al’pert [18 ], Fock [19], 



and others. The results of their analyses indicated needs for further refinement 
of the model such as the introduction of terrestrial sphere stratification, anistropy, 
continuous ionospheric stratification, polarization coupling and anisotropy at the 
ionosphere and the effect of ions, both heavy and light. 


2. Recent Investigations 

Among the most recent and active investigators are Budden [20], Wait [ 21] , 
Johler [22 1, and Kranushkin [23]. 

Wait is a proponent of the mode theory (or the residue theory). To account 
for the behaviors of the fields in the presence of an ionosphere, he considered 
successively the ionosphere as sharply bounded, gradual, and then stratified; 
isotropic and the anisotropic; uniform and non-uniform. The stumbling block of 
the mode theory was the transcendental model equation in Hankel functions of 
complex order for determining the modes of the propagating waves . In its solu- 
tions, various approximations were employed. In view of this difficulty, Johler 
attempted the solution of the VLF propagation problem using only the original 
zonal harmonic series. With the improved summation method, the series could 
readily be summed on a large computer in spite of a large number of terms 
needed. Budden, on the other hand, advocated a direct and numerical solution 
of the full-wave equation taking into account of the coupling given arisen by the 
anisotropic and non-uniform nature of the ionosphere. 

Kranushkin resolved the problem using the theory of non-self-adjoint dif- 
ferential operator eigenfunction expansion. He called his method the normal 
wave solution. He considered the ionosphere as stratified media with the bound- 
aries being weakly angle dependent. 

We will consider in the next section the four mathematical models currently 
being used; namely, 

i. the full -wave solution, 

ii. the mode theory, 

iii. the zonal harmonic series, 

iv. the normal wave method. 


3. Theoretical Background 


Most of the theoretical considerations of problems involving the propagation 
of radio waves are based on the assumption that the electromagnetic field satis- 
fies Maxwell’s equation in the form 
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V X E = - i a)fx Q H 
Vx H = icoeE 


In the above equations, a time harmonic factor exp(i« t) is understood; E and H 
are the complex representations of the electric and magnetic field vectors; /j-q 
and e are the permeability and permittivity of the material medium. For most 
non-ferromagnetic media, the permeability equals to that of the free space. For 
the ionosphere, the permittivity is a tensor. If e is put equal to e Q k, where e 0 
is the permittivity of the free space, then k , the dielectric tensor can be written 
as 



/ U 2 -^ 2 Y 2 

x - in YU - ImY 2 
\ im YU - InY 2 



U (U 2 -Y 2 ) 


in YU - ImY 2 

U 2 - m 2 Y 2 
- i4 YU - mn Y 2 


- im YU - InY 2 

i^YU - mn Y 2 
U 2 - n 2 Y 2 



The standard URSI notation is used here. X is the square of the ratio of plasma 
frequency to wave frequency; U = 1 - iZ, where Z is the ratio of the electron 
collision frequency to the wave frequency; Y (- 1 , m, n) are the components of the 
vector Y referred to some rectangular cartesian coordinates. Y is the magni- 
tude of the ratio of the gyromagnetic frequency of electrons to the wave fre- 
quency in the direction of the earth’s magnetic field. 

The dielectric tensor is asymmetric. It implies that the ionosphere is an 
anisotropic medium. The changes in the electron density and collision frequency 
and their variations with height contributing to the inhomogeneity to the property 
of the ionosphere of the VLF wave propagation. 

In all the problems considered in the sequel, the primary source is a Hertzian 
vertical dipole. 
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II. THE DIFFRACTION MODEL 


We will consider first the diffraction problem together with Van der Pol and 
Bremmer. 


An electric dipole is situated at Q , a distance b away from the center of the 
sphere of radius a . The observation point is at P, a distance r away from the 
center of the sphere. The points Q and P are separated by a distance R. 

The electromagnetic fields, E and H satisfy Maxwell's equations. A time 
factor exp (- i cot) is used. The analytic problem consists of finding an electric 
Hertzian vector potential function ttt, where r is a position vector, such that 
it satisfies the following differential equations and boundary conditions: 

(1) (V 2 + k 2 ) 77 = 0 (r > a), 

(V 2 + k 2 ) 77 = 0 (r < a) . 


( 2 ) 


"A 

_ __(r77) and k 2 77 are continuous at r = a, 
3r 


( 3 ) 


A singularity at Q of the form 
field. 


ikjR 

e 

ik 1 R 


which is further called the primary 


The electric and magnetic fields are given by 


E = k 2 


d 2 

dr 2 


H. 


0 r 3r3i9 


(r77> , 


^ CO GO 


and 


k 2 - 

K l,2 


"1,2 + iT l,2 " 
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1. The Zonal Harmonic Series 


The solution of the total Hertzian potential function in the form of the zonal 
harmonic series for r > a is, 


77 


ik.R 
_ e 

tot ~ lk"R 


T" 1 (2n + 1) R (k ‘ a) C<‘> (k,b) (k,r) P (cos 0 ) (1.1a) 

4* ^>(k,a) 


and for r < a is 


7T t 


tot 


k? 


n =:0 


(2n + 1) ( R n + 1)-T— T ^n X) (k 1 b)0 n (k 2 r)P n (cos 5) 

' n V K 2 3 ^ /I 


(1.1b) 


where 





^ 2) (x) = 


77 h< 2) 

2x n + 1/2 


(x) = x n 




'/'n (X) (X) + ^ 2) (X)> 


J , ,, (x) = x n 
2 x J n + 1/2 v ' 




(H = Hankel function, J = Bessel function) and where 
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i_d_ 
x dx 


log (x 4 n (x) } 


R =. 


-* x= kj a 


log {x 4 (x)} 

x dx 


-I x= k- 


L— 

x dx 


l°g {x C< :1) (x)} 


-» x= k, 


1 d 
x dx 


log {x4 n (x)} 


-J x=k 


It is well known that 


ikjR 


, kiR 


U V 

= (2n + l) .y 1) (k 1 r)T/r i (k 1 b) .P n (cos^) (r > b) 

n=0 


ikjR 00 


k r R 


= 2] (2n + l)£i^b)^r)P a (co.«) C* < b ) 

n=0 



By substituting (1.3b) into (1.1a), (1.1a) can be expressed explicitly as 
harmonic series. Furthermore, on the surface of the earth where r - 
be reduced to the following 


where 


77* 


tot 


(ka) 3 


L 

n = 0 


(2n + 1) 


O', a) 



(cos 8) 


N (x , y) = 


II. log (x (x)} 
x dx 


— * 


x= kj a 



In the reduction to Eq. (1.4), use is made of the Wronskian: 


( X > ^n X ) 00 - ^n X) < X > ( X > = ~T 

X 2 



— ( 1 . 2 ) 

2 a 


(1.3a) 


(1.3b) 


a zonal 
s , it can 


(1.4) 


(1.5) 


( 1 . 6 ) 
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2. The Watson Transformation and the Residue Series 


Following Watson, we transform (1.4) into an integral over n leading to 


77 


tot 


1 

(k,a) 3 



ndn t&/, ( k ,b) 1 „ 

cos' (mr) C (l> 2 (kia) N „_ 1/2 2 


[cos (77 — (9)] 


(1.7a) 


where Cj is the contour which enclosed the positive real axis in a counter- 
clockwise direction. By a judicial manipulation of C x , Eq. (1.7a) is equivalent 
to the following: 


77 


tot 


I 


ndn ^nli/2 1 


ki a)3 J L C aj /? ( kl a) N n-l/2 


n-1/2 


[COS (77- 0 )] 



f (n) dn 


(1.7b) 


where f (n) is the integrand of (1.7a). The line integral in (1.7b) is usually con- 
sidered negligible. L is the contour encloses all the poles in the integrand of 
(1.7a) other than the ones due to cos (n77). 

The first term in (1.7b) enables us to obtain a physical interpretation. Let 
us expand l/cos (n77) as 


1 

COS n77 


2e i7rn 



e i7Tm (2n + 1 ) 


( 1 . 8 ) 


We can write (1.7b) over the integration path L as 



(1.9) 
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where 


77 

m 


2e i7rm 


(k ia ) 3 


1 


ndn e i7rn ( 2m + 1 ) 


^Ii/2 »|b) 1 

'ill a N '"‘ 2 


x P n-l/2 t cos (rr-0)] 


( 1 - 10 ) 


The integral in (1.10) can now be represented as the sum of the residues of 
all the enclosed poles in L, and it is written as 


77 = 

m 


4 77 


i e i7rm V 1 

— r / ns< 

( k i a ) 4=o 




(2m + 1 ) ^j-1/2 (k l b) 

eiw, < k l a ) 


3N 


n-l/2 


P n -1/2 [COS (9 — 77 ) ] (1-11) 


in 


Upon using a well-known asymptotic formula for the spherical harmonics 
of complex order with positive imaginary part, (if 6 is not too near 0 or n), 


P , „ Tcos (n-M = 1 (1-12) 

n - I/a ✓2rniiin« 


^ 2/2ff e i7r ( m + 3/4) 

(k ia ) 3 •/ sin (9 


' L 

s = 0 


«;>,/, ( k .» 


n 


1/2 


BN 


n-l/2 


Bn 



in (# + 27Tm) 

e s 


(1.13) 
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The last factor in (1.13) becomes 

in ( 6 + 2-rrm ) ( i® . - /3. ) ( ® + 2nm) 

e 5 = e s 


(1.14) 


where 


n s = a s + 1 A 0 3 s >0 ) 

Hence (1.13) represents a wave travelling round the sphere m times before 
reaching the receiver with a phase velocity and attenuation determined by a s 
and fi a . Because of the fact that a wave travelling m times round the earth be- 
comes greatly attenuated, only the m = 0 term contribute significantly to the 
fields. Thus, 


^tot = 



A, Qlv* ( k . b > 

s=0 ^nl-1/2 ( k l a ) 


n 


1/2 


'3N 


(1.15) 


n s -l/2 


3n 


3. The Eigenfunction Method 

Van der Pol and Bremmer also consider the possibility of casting the radio 
wave propagation as an eigenvalue problem. In this case Eq.’s (1.1a) and (1.1b) 
are determined by the following considerations: 

(1) (V 2 + k 2 \ 2 )0. = 0 (r > a) 

(V 2 + k 2 X 2 . ) 4 >. = 0 (r < a) 

c) 

(2) k 2 <p. and — (r 4>.) are continuous at r = a . 

3 or 3 
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(3) k? 


i f 4 >)ar^l f* 

*'r=n «^r = 0 


4>? dr = 1 


where 0. ’ s are the eigenfunctions corresponding to the eigenvalues X.’s. The 
primary field pertaining to a point source leads to a three dimensional Dirac 

function at Q. 

k i 

The eigenfunctions are of the form 

'A nJ £ (1 > (X n>j k 1 r)P n (cos &) (r > a) 


r 


4> . = 
n. J 


B n , j <^n,j k 2 f ) P n ( COS ^ ( r < a ) 


(1.16) 


where n is a positive integer. The total field satisfying conditions (1) through 
(3) becomes 


CO n 


47ri \ 1 \ ' 

77401 " ~ Z_. z_. 

*• a : n 


—7/ n-k* .) 

n=0 j =0 v n , j y 


A "’ j « 15 <A n .,' k ib) 


£< 1} (X n j k 1 r)P n (cos 0) (r > a) 


(1.17a) 


00 n 


. .. iii y y jlli Bnj s<‘> (\ . k.b) 

"" k, L L, (1 _*a ) n 1 

1 n=0 j=0 n.j^ 


<£ n ( X n,j k 2 r ) P n < CO S ^ ( r < 8 ) 


(1.17b) 
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where 


A 2 . = • 
n. 3 


( 2 n + 1 ) 


27 rk 2 a 3 {£(*) (x )} 2 


(£j4st) l06 {xi » ,) ( x)> 


and 


B 2 . =- 

n. J 


_i /_!_ + iiL\ log {y^ n (y)> 

k 2 \dy 2 ydy/ 


( 2 n + 1 ) 


277 k 2 a 3 {yp n (y )} 2 

/ d 2 Id 


k 2 


+ l d\ 

_>? 

\dx 2 

X dx ] 


+ — t- ) lo g {y’/'n (y)> 
dy 2 ydy; 


(1.18a) 


(1.18b) 


The difference between the earlier method and the method of the eigenfunc- 
tions is clear. In the latter the summation is done over the integers whereas in 
the former it is done over the complex roots of the modal equation. 

9 

To sum up, we have seen in this section that there are now clearly three 
different methods used in the solution of the VLF propagation problems. They 
are (a) the zonal harmonic series method; (b) the residue series method; and 
(c) the eigenfunction method. In the next two decades, most of the analyses are 
done following these three methods. 


THE EARTH-IONOSPHERE MODEL 

To account for the effects on the VLF wave propagation due to the presence 
of the ionosphere, additional boundary conditions were introduced. The ionosphere 
was first considered as a sharply bounded homogeneous medium, then as stratified 
media. Magneto-ionic theory was also introduced to take care of both the iso- 
tropic case where the earth magnetic field was ignored and the anisotropic case 
where the earth magnetic field is included. The additional conditions introduced 
led to a great deal of complexity; nevertheless, the mathematical models em- 
ployed had their genesis in the diffraction model. 
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The Mode Theory 


Among the most active proponent of the mode theory is J. R. Wait. We shall 
look at the problem with him. 

A. Isotropic, uniform, sharply bounded ionosphere. The ionosphere is treated 
as a concentric homogeneous sharply bounded ionosphere. The source is a ver- 
ticle dipole. The fields are expressed in terms of a Hertzian vector which has 
only a radial component U . In the region between the surface of the earth ( r = a ) 
and the lower edge of the ionosphere (r = c, or 8 +h ), 


E 0 




(ru), 


1 B 2 


r Br ~dd 


(ru), 


- 1 6 0 ) 


>u 

Te 


H r = H, = 0. 


( 1 . 1 ) 


The time factor exp (icot) is understood, e and /jl are the electric constants 
and k = (e jj ) 1/2 co . Subscripts g and i are added to the field quantities where 
reference is made to the region r < a and r > c , respectively. To account 
for the singularity introduced by the source at r = b , the Hertzian potential 
satisfies the inhomogeneous wave equation, 


(V 2 +k 2 ) V = (1-2) 

2 tt r 2 sin 6 

for a < r < a + h , where the § 's are the dirac delta functions. The field in the 
region a < r < a + h is written as the sum of two parts U e and U s , where the 
latter is a solution of the homogeneous wave equation and the former satisfies 
the singularity due to the dipole source. 


^ (2q + l) CA q h< 2 > (kr) +B q j q (kr)] x P q (cos 0) (1*3) 

qssO 
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where j q (kr) and h^ 2 ) (kr ) are spherical Bessel function of the first kind and 
Hankel function of the second kind respectively; and P q (cos#) is the Legendre 
function. And 


r 


u 

C 47 T 


(2q + 1) P q (cos 6) J 

q^o 


j (kr) h (2) (kb) for r < b, 


h< 2 > (kr) j q (kb) for r > b, (1.4) 


The Hertzian potentials U g and U. are solutions of the homogeneous wave 


equations 


(V 2 + k 2 ) U g = 0 for 0 < r < a 


and 


(1.5) 


(V 2 + k 2 ) U i = 0 . for r - c. 


To satisfy the boundedness of U at r = 0 and r = », the solutions are 

= £ (2qtl)P < (C0S *>■« <1-6) 
q 

and 

U i ^](2q + i) P q (cos 9 ) b q h< 2 > ( k i r ) (1-7) 

q 

where a and b are constants of integration. For the boundary conditions, 
Wait introduced the impedance concepts; thus, for q f th term in the series they 
are 

E^ q) = -Z< q > H£ q) at r = a ( 1<8 ) 

E^ q) = Z( q) H^, q) at r = c ^ 1 ’ 9 ^ 
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where 


and 


Z (q) 

g 


i ^ [riq(kq ° ] 
ieu r J q ( k q r ) 


JL [rbf)( V )] 

2<q) = - 1 dr 


1 eco 


r b< 2 > (k.r) 


These can be rewritten upon replacing kr by x as 


1 JL (XU) = i (Z< q )/7?)U for x = ka 
x 3x g 


-L -L (XU) = - i (Z < q Vt?) U for X = kc 

x Bx 1 


The solution for a < r < c is 


where 



CO p 

y 1 ' (2q + 1) h< 2) (kb) hj 1 ) (kr) x P q (cos 9)-l 

q=0 


1 + R< q > 


hq 1) (ka)h< 2 > (kr) 
h< 2) (ka) (kr) 


D =l-R< q )R< q > 

q g 1 


h (2) (kc)h< i:) (kb) 

l 1 p(q) q q 1 

1 h< 1J (kc)h< 2) (kb) 


(ka) h^ 2 ^ (kc) 


h< 2) (ka)h^ 1 ) (kc) 


R( q ) = 


_d_ 

dx 

IT 

dx 


«u x hj 1 ) (x)] x=ka 
Wu x h< 2) (x)] x=ka 


- iZj q )/7] 

- iZ^/r] 


( 1 . 10 ) 

( 1 . 11 ) 

( 1 . 12 ) 

(1.13) 

(1.14) 

(1.15) 

(1.16) 

(1.17) 
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(1.18) 


g- h $ 2) < x ) ] x =kc + iZi q) /T7 

-T- t^ u x h J l) (x)lxpkc + iZ i q) /T? 


By means of the Watson transformation, the solution becomes 


v — i v + 2 F 

U = - ike / — — — — (kb) h^ 1 ) (kr) — P [cos (w - 0)] (1.19) 

L i sin vv n' 


where = B/9v D v . All the values in the summand are to be evaluated at the 
poles of f ( v) which are the roots of the model equation 

D v = 0 (1.20) 


If use is made of the Airy integral as expounded by V. A. Fock in place of the 
Hankel function, the radial component of the electric field can be written in the 
form 


E r £ 


. - ika& 


a (<9 sin (9) 1/2 


(1.21) 


apart from a constant factor. 
And 


V = e 


iw/4 


n' /2 <£ e '“ t [w ' (t - y)+A(t)w ^ (t - y)1 dt (i.22) 

^ 77 ' 7 [ w ; (t) - qwj (t) 3 [1 - A (t) B (t) ] 


where 


x = (ka/2) 1/2 6, y = (2/ka) 1/3 k (r - a) 


(1.23) 


Now the contour is to enclose the complex poles t n’s which are solutions of 

1 - A (t) B (t) = 0 (1*24) 
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where 


A (t) = - 


w 


[ (t-y 0 ) + qjWj (t -y 0 ) 


w 2 (* -yo) + qi w 2 C* -y 0 ) 


(1.25) 


B (t) = - 


w' (t) - qw 2 (t) 
w| (t) - qwj (t) 


(1.26) 


y 0 = (2/ka) 1/3 kh, q. = -i (ka/2) 1/3 Z/t) 0 , t] q = 120 tt 


and 


q = - i (k a/ 2) 1/3 


ie o" 


cr +16 co 
g g 


1 -• 


1 € o a> 


, 1/2 


a g +1£ g W 


The residue or mode series is 


,/o e‘ ixtn [wjCtn-y) +A(t n ) w 2 (t n -y>] 

V = - 2 (-77 x) 1/2 e* 1 /4 y — i (1-27) 

ns 0 K (t n ) -qw x (t n )] 


dt 


A (t) B (t) 


t=t 


It may also be written as 


v - 4 ^*) 17 - 2 e^ 4 Y " e' iXt "G n (y)G n 

y o 4^0 


(y) An 


(1.28) 


where y = (2/ka) 1/3 kZ 0 . The functions G n are height functions and they are 
normalized to unity for y or y 0 equal to zero. Explicitly 


c n (y) = 


f (yy) 


(1.29) 


f (t n , y) = w x (t n - y) + A (t n ) w 2 (t • - y) 


(1.30) 
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The function A n is a modal excitation factor. It is normalized so that, in the 
limit of zero curvature and perfectly reflecting boundaries, it becomes unity for 
all modes. Explicitly, 



^n-yp-qf) ^2 (<„) -q w 2 0^)3 2 
Cw 2 (V-y 0 ) +q i W 2 ( t n-yp ) ]2 


(1.31) 


b. Anisotropic Ionosphere. For the case of an anisotropic ionosphere where 
the earth magnetic field is not ignored, the electromagnetic field components 
in the homogeneous space (a - r 5 c) are expressed in terms of purely TM and 
TE waves. They are derived from two scalar functions U and V as follows: 


E ' = (|? +k2 ) <rU) 


E* = 


1 B 2 


(rU) - 


i /xo> B 


r B^Br r sin 9 B<£ 


(rV) 


1 B 2 
^ r sin 6 B0B r 


Ej. = - A - . . (rU) + JL (rV) 

r do 


(1.32) 


H ' = {£ + k2 ) (rV) 




1 B 2 


(rV) + 


i e co B 


r B(9 Br r sin 9 B 


(rU) 


% = Ar ( fV ) - — ( rU ) 

v r sin 6 or ocp r do 


The functions U and V satisfy 


(V* + k 2 )|“j = 0 


(1.33) 
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i 

in source free region. The source again is a vertical dipole. The primary field 
of such a source is of a purely TM character. It is derived from a single scalar 
function U e satisfying the following inhomogeneous wave equation, 

(V 2 + k 2 )U e = c 3(r - b > 3(g) (1.34) 

2 tt r 2 sin 9 

The boundary conditions, in terms of surface impedance, are 


|-(xU) = iA g xU 


JL(xV) = (i/A g ) xV 


x=ka 


(1.35) 


where 


A g = Z g /T 7 , U = U e+ U, V = 17 V; 




and 


i-(xU) =A U ^(xV) - iA 12 (xU) 


“i (*V) =A 21 2(xV) - iA 22 (xU) 


:X=kc 


(1.36) 


where 


and 


Aj^ — Z^/t], Aj 2 - Z 12 /rj, A 21 Z 21 ^^ 


^22 _ ^22 // ’^ - 


The radial field in the concentric waveguide region written in the matrix 
form is 


" E r _ 

~ e -ika ^ 

'F e “ 

i [L Q Cl) I d S 

_^H r _ 

a (# sin #) 1/2 

_ F h_ 

477 


(1.37) 
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where 


L F hJ 


= -2(ttx) 1/2 e“ i7T/4 — 


-ixt n [w,(t n -y) +A(t n )w 2 (t n -y)] 


[w'(t n )-qw 1 (t n )3 


Bt 


A(t) B(t) 


(1.38) 


-J t 


The roots t n f s are determined from the matrix equation 

[A (t ) ] [B (t ) ] = I 


(1.39) 


where 


and 


A (t) = 


,Ri 


i R n 


n R i i R i 


exp [-1 (4/3) (y 0 - 1) 3/2 - in/ 2] 


B (t) = 


w' (t) - qw 2 (t) 
wj (t) - qwj (t) 


w 2 (t) - qw 2 (t) 
(t) - qwj (t) 


(1.40) 


(1.41) 


q - - i (ka/2) 1/3 A g , with A g = 


t ie o" N 

1 W2 I 

\°g + i V"/ 

1 \ 


1 - 


ie 0 a> 


an 


cr + i e co 
g g 


and 


q = - i (ka/2) 1/2 A h , with A^ = ( _ 1 £ 1 

8 8 ' ie Q a> 


1/2 


ii r ii ’ i R n» ii r i > and jR^ are the reflection coefficients , where II and 
refer to the direction of the electric vector, relative to the plane of incidence, 
of the incident (first suffix) and reflected (second suffix) wave. 
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The Zonal Harmonic Series 


The stumbling block of the mode theory is the rather formidable transcendental 
modal equation to solve. Recently Wait and Spies made detailed computation of 
the characteristics of the least attenuated modes. However, a single mode is not 
adequate to describe the complete field satisfying Maxwell's equations at either 
intermediate and great distances for daytime models of the lower ionosphere. 

In the case of a nighttime model even more modes are needed at all distances. 

With the speed of the modern digital computers, the argument once held 
against the zonal harmonic solution of the YLF propagation problem is no longer 
convincing. Granted that the zonal harmonic series converges rather slow. 

Johler showed that with improved summation technique, the results of the compu- 
tations indicate that full rigor can be achieved with comparative ease at fre- 
quencies less than about 50 kc/s, leaving only the assumed model for the trans- 
mitter and the propagation medium and avoiding the complications of the Watson 
transformation. 

For a simple model of the terrestrial sphere of radius a surrounded by con- 
centric ionosphere from r = c to r = ®, the field components are expressed in 
terms of the Hertzian potential functions, v e and n m as follows, 

E = — — 

r rb sin 9 


E e = — — 

9 rb Br 3(9 


3_ 

30 


sin 0 


377 e 

30 


( f <TT ® A 


_ 1 k 2 37T e 
^ b /j- 0 i co 35 


v _ V-o 1 3? r m 

* b ~ 


( 2 . 1 ) 


„ - 1 'd f .. a d 7T m 

H — — - sin 0 

r rb sin 0 30 \ 30 
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where tt 9 and rr m satisfy the wave equations 


(V 2 + k 2 ) 


TT e 
7 T m 


=: 0 


( 2 . 2 ) 


The media of propagation are characterized by their electric constants in 
the wave number, k . Thus, for air, with index of refraction, 77 , or dielectric 
constant e 1 = 77 2 


or 


k 2 = k 2 =f_e t 

e 2 



For the ground 



(2.3) 


(2.4) 


where cr is the ground conductivity and e 2 is the relative dielectric constant. 
The wave number of the ionosphere is rather intricate. 


k 



(2.5) 


There are four distinct propagation components with the complex index of re- 
fraction. 77^ corresponding to ordinary wave (0), and extraordinary wave (e), 
upgoing wave (i), and downgoing wave (r). They are related to the roots of a 
Booker type quartic in the parameters, £ , which is the result of the simultaneous 
solution of the Langevin equations of charged particles and Maxwell's equations. 
The index of refraction, 77 , 


77 2 = £ 2 + sin 2 0 

where <f> i is the angle of incidence on the lower ionosphere plasma. 


( 2 . 6 ) 
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The boundary conditions are 


k l Oo + 77 1> = k 2 OP 



T77P 


(rr/p 


77 ™ = 77 ™ 


(2.7) 




r - 


a 


k l K + W I) = ^0 1 " T 37- (r 


r or 


77: ( r7T S + r7T i) - ^0 1 w Q^e ^3 


77™ = Q — JL (r 77®) 
1 ™ e r dr 3 ' 


( 2 . 8 ) 


^0 iM ^( r7T l) = Qem k 3 77 


where in the quasi-longitudinal approximation, 


k 3 /k l 

^5 /kg/kj - ~p 


E - sin <£. 



(2.9) 


/? = sin 2 4> 2 is the complex direction sine squared of propagation, and the minus 
sign (-) is taken with the ordinary wave, k 3 0 and the plus sign (+) is taken with 
the extraordinary wave k 3 e . The solutions are 
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J2 (2n + 1) C n (kjb) yp n (k t r) P n (cos (?) , 

n=0 


(r <b) 


1 

rb 


( 2n + 1) l n (kj r) ^ (kjb) P n (cos 0) , 

naO 


(r > b) 


V \ = F7b £ (2n + 1} [b " ^ < k i r ) + C n '/'n ( k X r >] P n (COS 5) 

1 n=0 


77 


m 

1 


l 

kj rb 


E 

n=0 


(2n + 1) [b” C n (k x r). + c” C n <k x r)] P n (cos 0) 


00 

— T” 1 (2n + 1) a® \p n (k 2 r) P n (cos 6) 

k 2 n=0 


77 ” = 


1 

k|7b 


w 

L 

n=0 


(2n + 1) 


a”>A n (k 2 r)P n (cos 6) 


77 


e 

3 


1 

kfTU 


00 



(2n + 1) d® 5 n (k 3 r) P n (cos 0) 


7T* 


1 

kfTb 


L 


(2n + 1) d” £ n (k 3 r) P n (cos 0 ) 


( 2 . 10 ) 


The constants a® , b„ and the like are determined by the boundary conditions. 
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For a stratified ionosphere model, the matching for boundary conditions 
could be a rather formidal task. (For details see Reference 22.) 


The Full-Wave Theory 


It is seen that when coupling is introduced into the wave propagation prob- 
lem, no matter what mathematical model one uses whether it is based on the 
residue series or the zonal harmonic series , one often resorts to the digital 
computer for numerical computations. Budden is one among many of the English 
schools who rely more and more on the direct solution of the full wave equation 
using computers . This approach is formally of complete generality and does 
not depend on the assumption of a special ionospheric model. 

In deriving the differential equations, Budden uses the cartesian coordinates 
with the z-axis directed vertically upwards. A plane wave is incident at angle 
<9. to the verticle on the ionosphere from below. Let s • sin 6> i , c = cos 9 { . 
Let the ionosphere be stratified such that in each of which the medium is con- 
sidered homogeneous. For the waves in each stratum B/Bx ( ) - - iks( ); 

B /By ( ) = 0. Then Maxwell’s equations give, 


dE X 

— + iksE = - ikH 
dz z y 


- iksE = - ikH 
y z 


(3.1a) 


dH : u 
l = _D 

dz 6 0 * 


dH 


iksH = — D 
dz e Q y 


- iksH = —B 

y e 2 
c o 


(3.1b) 
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From the magneto-ionic theory, after the Z -components of the fields are 
eliminated, we can write the following matrix equation, 


where 


de 

dZ 


- ikTe 


(3-2) 



(3-3) 


and T is a 4x4 matrix whose elements are related to the parameters describing 
the ionosphere, such as the electron density, the electron collision frequency, 
the earth’s magnetic field and others. 


The matrix T has four characteristic roots or eigenvalues q. ( i - 1, 2, 3, 4) 
which satisfy the characteristic equation 

det (T - ql) = 0 ( 3 - 4 ) 


where I is the 4x4 unit matrix. 


In order to solve the matrix equation, it is convenient to introduce a new 
column matrix f , thus 



(3.5) 


such that 


e = S f (3.6) 

where S is some 4x4 matrix whose elements are functions of z to be determined. 
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Assuming that S is non-singular so that its inverse S _1 exists. Substitution 
of (3.6) in (3.2) and premultiplication of S _1 gives 

f ' + ikS' 1 TSf = - S" 1 S' f (3.7) 


where the dash indicates the derivative with respect to z . 

The S is so chosen as to make f . the only element of f appearing on the 
left hand side of the i’th equation in (3.7). This means that S' 1 TS must be a 
diagonal matrix, this can be done if the characteristic roots q t of T are all 
distinct. It can then be shown that 


S- 1 TS = 


0 

0 

0 


0 

q 2 

0 

0 


0 

0 

% 

0 


0 

0 

0 

% 


= A 


(3.8) 


and Eq. (3.7) becomes 

f' + ik Af = -S- 1 S' f (3-9) 

Equations (3.9) are the four, first order, coupled differential equations. For 
a homogeneous medium S' = 0 and (3.9) becomes 

f ' + i k A f = 0 (3.10) 


or 


f ( + ikq. f . = 0 (i = 1 , 2 , 3 , 4) ( 3 * n ) 

The independent characteristic waves in a homogeneous medium are determined 

by 


f i =e~ lkqi * (i = 1 , 2 , 3 , 4) 


(3.12) 


27 



Now, premultiply ( 3 . 9 ) by any 4x4 matrix, say F, (3.9) then can be written 
as 

(Ff)' + (ikFA-F')f = -FS -1 S' f (3-13) 

Choose F so that 


F' “ ikFA, 


(3.14) 


and 


(Ff)' = -FS -1 S' f 


(3.15) 


and integrate (3.15) in the form 


f z 

f = F -1 - F -1 FS _1 S' f dz (3.16) 


From (3.14) F may be taken as the diagonal matrix whose diagonal elements are 
exp (ik J z dz); hence, if f 0 is the fundamental solution of (3.10) namely 



(3.16) becomes 


f = f. 


f" 1 
1 o 


S'f 0 dz 


(3.17) 


(3.18) 


The evaluation of S _1 S' and its properties can be found in Budden’s book. 

The thing of interest to us here is that we have a set of first order differential 
equations (3.15), which can be integrated numerically choosing various iono- 
spheric models. A work of this kind has been undertaken by Budden and his co- 
workers in the Cavandish Laboratory. Their method of attack is to assume some 
plausible law for the variation of electron density and collision frequency with 
height in the ionosphere, and to calculate the reflecting properties of this model 
ionosphere for long and very long waves. This is repeated for a series of models, 
and those models are selected which have properties most closely resembling 
the experimental results. 
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The Normal Wave Solution 


In all the previous methods of solution of the VLF wave propagation problem, 
numerical checks were made by assuming a prior knowledge of the parameters 
used. The selections of the parameters though based on sound and scientific 
judgment, were, to say the least, arbitrary. Furthermore, these methods could 
account for only homogeneous tracks (night or day) up to now. The diurnal 
changes in both the phase and the amplitude of the VLF waves remained mostly 
unexplained. In view of the inadequacy in these approaches, Krasnushkin proposed 
a mixed method. He suggested to use the data known in the short range propaga- 
tion to evaluate some of the unknown parameters. This is basically an inverse 
method. With the unknown parameters determined, a direct method is used to 
evaluate the fields at great distances. In doing so, some of the difficulties en- 
countered in evaluating the roots of a transcendental equation as in the mode 
theory could be circumvented. 

The mathematical models he used are based on the theory of non-self adjoint 
differential operator eigenfunction expansions. He called these functions normal 
waves. 


a. Stratified ionosphere. In a spherically stratified media, let the dielectric 
tensors be as follows, 



e k 

0 

0 


r r 



11 

__ 

ID 

0 

e k 

^ee 

F k 

e e<p 


0 

e k 

e 66 


(4.1) 


r k _ 1 £ r < r R k = 0, 1 , 2 • , N 

where, for the zero layer (k = 0), or (r x = 0, r Q ) 

4 77CT n 


-0 _ ,0 _ ,* , , 
rr ~ e ee - 6 + 1 


0 i eL = o 


CO 


z e<t> 


(4.2) 


for the first layer (k = 1); (r 0 = a, i^) 


€ rr = e ee - 1 e d<t> ~ 0 


(4.3) 
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and for the other n - 1 layers (k = 2, 3, . . . . N) 


( r i> r 2 ) ; < r 2> r a) 

the tensor ||e || is as follows, 


( r N-l ’ r N _ °°) 


'■ _ 1 


e ee - 1 + ■ 


iaJ 0 ( V e ff “ ia) ) 


CO t(v eff - i oS ) 2 + 0 ) 2 ] 


(4.4) 


e e<t> ~ 


iaJ 0 W H 


CO [<v eff - i CO) 2 + co 2 ] 


where 


0)2 = 


47rNe 2 


m 


eH„ 


0> H - 


me 


(4.5) 


The components of the field in each layer are 


1 3 / • _ i-CO u u> k , ^77 T k 


r sin 6 dt 9 


(sineHj) = -— e r r 1^ , 


7 37 <r “ T 1 [e «*S + 4 * E P. 
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— ^ e \e E 5 + e 44> ’ 
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r sin 9 B 9 


k\ _ la> jjk 


(sin 0 E{) = - ^ 


_1_B_ 

r 3r 
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1 

r 



BE^ 

~W 



(4.6) 


They are subject to the following conditions: 

1. Tangential components E 0 , E^, H 0 , are continuous on the spherical 
surfaces r R = const., R = 1, 2, . . . , N. 

2. In the center of the earth (r = r =0) the fields are bounded and 

3. At infinity they tend to zero. 

Now, introduce the Hertzian potential function 
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where k 0 = co/c is the wave number of the free space. 
Eq. (4.6) can be expressed in terms of A and B 
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The continuity conditions for E g , E d , E^, at r = r k become 
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Now decompose (4.7) into the product 


B 


Y (r ) 

A 


Z (r) 


'P (0) 


(4.14) 


Then will act only on I Y I and l e only on 6. After replacing B/Br and 
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d/dr one obtains the differential operator L f defined by (k). Similarly, one 
obtains L 0 defined by ^ e (k). Then one has a single operational equation: 
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To find the solution of (4.15), we will find first the solution of the homogeneous 
equation 
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From (4.14), we have 
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L i9 'A (0) + x \p (6) = 0 
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The solution of (4.17 
ing to the eigenfunctions 
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Equation (4.18) is the usual Legendre equation 
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For the case of a verticle dipole source placed ar r = bwe may write the 
solution of (4.15) as follows, 
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where P is the dipole moment of the source and 
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Using the asymptotical approximation of the Legendre function, and expanding 
1/sin v. 7T , we write (4.21) as follows, 
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The separate terms in the braces of (4.23) denote the signals reaching the re- 
ceiver after passing round the earth n times. Because of the attenuation only 
the waves with n - 0 arriving at the receiver along the shortest path should be 
considered. Thus, 
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The field components are 
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For calculating the field in the atmospheric layer, a further simplification 
of formula (4.2) is needed. For the details involved, see reference [23]. 

b. Stratified medium with relief . To account for the diurnal variations of 
the field, Krasnushkin proposed another model to take into consideration of the 
slow variations of the dielectric tensor (4.1) along the path of radio waves. He 
assumed that 
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where 90° - 6 and qb are latitude and longitude and (j. is some small parameter. 
The operational equation now takes the following form 
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This equation is not separable due to dependence of L r on 6 and <fi. However, 
for small /x one may expand L r in powers of /. x. The first approximation yields 
the following eigenvalue equation: 
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From whence one obtains the set of eigenvalues x. (/x6, fifi), j = 1, 2, ... 00 
corresponding to the eigenfunctions 
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For <F (dfi) one obtains the equation 
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If one seeks the solution in the form 
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where U. is a slowly varying function of 6 and <fi. After substituting (4.32) into 
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(4.31), one finally obtains for the following expressions: 
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where normalization factors N. (0) and N. (P ) are evaluated at the sender and the 
receiver points. The path of the wave is defined by the eikonal function s , evalu- 
ated from the equation 
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If the geodetic line is followed, then 
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After substituting (4.36) into (4.33), one finally obtains the radial component 
of the electric field 
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SUMMARY 


So far we have seen four mathematical models in the study of the VLF wave 
propagation; namely (1) waveguide mode theory, (2) zonal harmonic series, (3) 
eigenfunction method, and (4) full wave theory. Of the four, the last stresses 
more on the computational aspects of the problem whereas the first three give 
analytic expressions as solutions of the problem. Models one and three are more 
elegant in form, each is capable of explaining the propagation problem in simple 
and tractable terms using mostly one or two modes. Model (1) is based on the 
concept of the classical residue series, whereas model (2) is formulated on the 
modem theory of the non-self adjoint operators. The stumbling block in the 
mode theory is the rather formidable transcendental equation whose solution 
is necessary to account for the various modes needed in evaluating the electro- 
magnetic fields. There is a corresponding transcendental equation in the eigen- 
value method, however, Frasnushkin showed that the difficulty can be circum- 
vented by means of what he called the mixed method; whereby he used the data 
of the short range fields and some of the known parameters to determine the 
unknown parameters . Then he used the known and the determined parameters 
to evaluate the long range propagation phenomenon. 

In contrast to the mode theory model and the eigenfunction model, the zonal 
series model relies only on the original harmonic series for its solution. The 
speed of the electronic computers offsets the objection of the direct summation 
of the slowly convergent harmonic series. A complete solution is thus made 
possible leaving only the assumed model for the transmitter and the propagation 
medium. The advantage of a descriptive discussion that a simple formula af- 
fords is lost. Perhaps when all the details of the attributes of the ionosphere 
are included, graphical solutions based on numerical computation are the only 
answers. In which case, the full wave theory model as is handled by Budden may 
be the best approach. Unfortunately, Budden* s is a planar model. To satisfy 
the long range propagation, the earth's curvature must be taken into consideration. 

The abundant experimental materials cumulated over the past fifty years in 
the field of VLF propagation have established many facts that the current theories 
are still incapable of explaining, such as 

1. the diurnal and the seasonal variations of the long distance VLF field; 

2. the dependence of both the field intensity and the phase on a heterogenous 
path of various illuminations; 

3. the effect of the solar flare on both the field strength and the phase. 
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It appears that some work could be done on the formulation of a new model 
which does not depend only harmonically on time in order to resolve some of 
the aforementioned observed facts. 
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A THEORETICAL STUDY OF ELECTRIC FIELDS 
IN THE LUNAR ATMOSPHERE 

by 

Richard A. Elco* 


ABSTRACT 

Lunar atmospheric electric fields are a result of charge sepa- 
ration and V x b interactions of the interplanetary magnetic field 
and the lunar ionosphere which is created by photoionization and 
charge-exchange ionization of the accreted neutral atmosphere. 
The interplanetary magnetic field provides a mechanism for trans- 
fer of momentum from solar wind protons to slow ions which in- 
creases the loss rate of accreted ions from the lunar ionosphere. 
A collisionless hydromagnetic model of the partially ionized gas 
flow on the moon yields the ordinary collisionless hydromagnetic 
fluid equations. Since the solar wind velocity upstream of the moon 
is supersonic, a shock wave forms on the moon. However, for a 
steady state shock, the shock must be aligned with the interplane- 
tary magnetic field. The possible geometry of the shockwave is 
discussed. The maximum increase in the total density and mag- 
netic field across the shock is three times the solar wind values. 
Therefore, the lunar atmosphere is extremely tenuous, having a 
total particle density of no greater than 50/cm 3 . Space charge 
electric fields in the lunar ionosphere are two orders of magnitude 
smaller than the V x B electric fields which are the order of ten 
millivolts per meter. 


INTRODUCTION 

The sources of the lunar atmospheric and surface electric fields can be 
grouped into the following three categories. The first, is charge separation in 
the lunar ionosphere. This mechanism is of primary importance if the iono- 
sphere is in diffusive equilibrium and the charge separation is due to the different 
scale heights for the ionic and electronic components of the ionosphere. The 
second source is a dynamic or V x B interaction which results from the motion 
of the lunar ionosphere in the interplanetary magnetic field. The third source 
is a surface electron sheath which exists to balance the solar wind current inci- 
dent in the lunar surface with the photoelectron current emitted from the lunar 
surface (Singer and Walker (1962)). The nature of the lunar electric fields is 
then directly related to the dynamics of the lunar ionosphere. The primary 
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physical processes which maintain the lunar atmosphere and ionosphere are 
(1) surface recombination of solar wind protons incident on the lunar surface 
and subsequent thermal emission of neutral hydrogen atoms, and (2) photoioniza- 
tion and charge exchange ionization of the accreted neutral hydrogen atoms. The 
lunar atmosphere is essentially a collisionless mixture of energetic solar wind 
protons, energetic electrons, thermal hydrogen ions, thermal or low energy 
electrons, and thermal hydrogen atoms. 

Model lunar atmospheres and ionospheres , in which the effects of the inter- 
planetary magnetic field in the solar wind are neglected, have been developed 
by Weil and Barasch (1963) and Hinton and Toeusch (1964). These models pre- 
dict neutral and ion particle densities several orders of magnitude greater than 
the solar wind density within one to two thousand kilometers of the lunar surface. 
The loss mechanisms considered in these models were thermal escape and 
elastic scattering by solar wind protons. However, as was pointed out by Harwit 
and Hoyle (1962) and Bernstein et. al., (1963), the interplanetary magnetic field 
can provide a mechanism for transferring momentum from the solar wind pro- 
tons to the slow ions generated by photoionization and charge exchange and hence 
increase the loss rate of ions from the lunar ionosphere. Other studies by 
Beard (1965) (1966) indicate that a collisionless hydromagnetic shock wave is 
formed by the interaction of the solar wind and the interplanetary magnetic field 
with cometary atmospheres which are, in some respects, not vastly different 
from the lunar atmosphere (Gold (1966)). Therefore, the presence of the 
interplanetary magnetic field in the solar wind can increase the rate at which 
ions are swept out of the lunar atmosphere and also give rise to the formation 
of a hydromagnetic shock wave which can drastically affect the lunar atmosphere. 

If the solar wind together with the lunar atmosphere is considered as a mix- 
ture of interacting particles, where the interaction is due to the ionization of the 
flow neutral hydrogen atoms and subsequent transfer of momentum from the 
solar wind protons to the ions by the magnetic field, then a collisionless hydro- 
magnetic model is appropriate for the description of the flow in the vicinity of the 
moon. In this analysis it is assumed that the moon does not have a significant 
dipole magnetic field. A multi -component hydromagnetic fluid model is de- 
veloped and the hydromagnetic jump conditions are presented. Since the solar 
wind velocity far upstream of the moon is supersonic, it is likely that a stationary 
shock wave does form on the moon. The consequences of the shock wave are 
discussed, especially as it affects the particle densities and electric fields in the 
lunar atmosphere. 


The Collisionless Hydromagnetic Fluid Equations 

The collisionless hydromagnetic fluid equations which are given below are 
derived in the same manner as for the ordinary two fluid (electrons and ions) 
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collisionless model except that, in this model, momentum exchange by processes 
other than collisions is taken into account. The exact process by which momen- 
tum is exchanged need not be specified, however, some possible processes are 
wave interactions between slow and energetic ions due to a two stream instability 
or magneto-ionic wave interactions due to local disturbances of the magnetic 
field when photo or charge exchange ionization occurs. Since the gas mixture 
is made up of energetic and slow electrons, ions and neutrals, a three component 
mixture can be used. The gas is assumed to be electrically neutral. 

The continuity equations of each component are - for the ions 
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The is the photoionization rate. No terms involving charge exchange ioniza- 
tion are involved because the net creation rate of neutrals or ions by charge ex- 
change is zero. The velocity w K is the mass center velocity of the K th com- 
ponent, i.e., n w = n s V s + n* V* where the s and asterisk refer respectively 
to the slow and energetic electrons and V is the particle velocity. A combination 
of Equations (1) through (3) yields the continuity equation for the fluid as a whole, 
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The equation of motion for the fluid is then 
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where P is the total pressure, i.e., 
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This form of the equation of motion which includes the neutrals is obtained for 
the following reason. During the charge exchange process a slow neutral atom 
exchanges an electron with an energetic proton which then becomes an energetic 
neutral, therefore the result is a net increase in momentum of the neutral com- 
ponent and a net decrease in momentum of the ion component (Biermann (1966)). 

The generalized Ohms Law for the fluid (Holt and Haskell (1965)) is 


^ t p o e /— ♦ — * 1 — ♦ — ♦ 

Lvp - _! [ E + u x B — J x B 

3 1 m e e m e \ n e e 



( 9 ) 


where J is the current density, e the magnitude of the electron charge and k i 
and ~K e are the rates of change of momentum of the ions and electrons due to 
interaction with neutrals . The rate of change of the ion momentum k i is nega- 
tive. It is also most likely that the electrons also experience a net loss of 
momentum during each charge exchange which is the order of m e / m. A £ . There- 
fore it is assumed that the R.H.S. of Equation (9) is zero which implies that the 
momentum exchanged in this way does not contribute to a current in the fluid. 


The energy equation for the fluid is 
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in which D/Dt = B/Bt + u • Vand e is the specified internal energy. If it is 
assumed, following Chandrasekhar (1960), that only the energy of charged particle 
motion transverse to the magnetic field contributes to the internal energy of the 
ionized portion of the fluid, then the specific internal energy e is 

n.WI+n 1 »i+|n,kT. 
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where W ± is the energy due to motion transverse to the magnetic field in a frame 
of reference stationary with respect to the mass center fluid velocity u, T a is 
the neutral gas temperature and k the Boltzmann constant. It is convenient to 
rewrite (11) as 


where 
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The partial pressures of the ion, electron and neutral gases are 

P. = n. W} , P = n Wj and P a = n kT 

i i X e ex a a a 

Then, from Equations (8) and (13) the total pressure is 
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The energy equation can then be written as 
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which indicates that the effective ratio of specific heats is somewhere between 
2 and 5/3 depending on the value of T a with respect to the values of nWj^ for 
the electrons and ions. 


The Lunar Hydromagnetic Shock Wave 


Studies of collisionless hydromagnetic shock waves have been made using 
nonlinear wave mechanisms and various forms of the hydromagnetic approxima- 
tions of the plasma equations (Beard (1965), Sen and Spero (1967)), to describe 
the shock structure and the initial and final plasma states. The hydromagnetic 
approach is used in the following analysis. While the hydromagnetic model does 
not provide information about the details of the shock wave structure, it does 
yield, in a straight forward way, the necessary conditions for shock wave for- 
mation and the jump equations across the shock wave. 
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Using the continuity equation, (4), the equation of motion (7)^the generalized 
Ohms Law (9) and the field equations V x E - - BB/Bt and V x H = J to de- 
termine allowed discontinuities in the fluid flow yields the jump conditions de- 
rived by Sen and Spero (1967) which are 
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where the subscript n refers to the component of a vector variable normal to 
the shock wave while t refers to those components which are in the plane of 
the shock wave. The symbol [ 3 is defined as 


[f (n)3 = f (0 + ) - f (OJ 

where the plane of the shock wave is located at n = 0. 

As a consequence of Equations (21) and (22), Sen and Spero (1967) show that 
no steady-state hydromagnetic shock wave can exist if B n 4 0. This implies that 
the surface of a steady-state hydromagnetic shock wave must be aligned with the 
magnetic field. Therefore, a shock wave formed by the flow of plasma around a 
sphere in a magnetic field perpendicular to the free stream velocity will be semi- 
cylindrical instead of hemispheric in shape where the axis of the semi -cylinder 
is aligned with the magnetic field vector. This conclusion agrees qualitatively 
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with the results obtained by Ludford and Singh (1963) in their study of the crossed 
field flow of a perfectly conducting, incompressible fluid over a non-conducting 
sphere. However, it is obvious that the semi -cylindrical shock surface cannot 
extend indefinitely on each side of the sphere, therefore, non-steady shock waves 
must exist on each side of the steady shock wave as is illustrated in Figure 1. 



Figure 1 -II lustration of combination steady state and non-steady state hydromag- 
netic shock wave on the moon due to interaction of lunar ionosphere with solar 
wind and interplanetary magnetic field. 


In view of the complexity and uncertainty implied by the above discussion 
of the shock wave geometry, the following assumptions are made in order to 
obtain an estimate of the post shock conditions which are relevant to the lunar 
electric fields. First, the interplanetary magnetic field is assumed to be per- 
pendicular to the free stream velocity of the solar wind and secondly, the steady- 
state shock wave surface is assumed to be a semi -cylinder aligned with the 
magnetic field as is shown in Figure 2. The local shock coordinates are the 
cartesian system x , y and z and the subscripts 1 and 2 refer to the pre-shock 
and post-shock states respectively. The magnetic field normal to the shock is 
then zero, i.e., B„ = 0, and also B = 0. 

x y 
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Figure 2-Simplified geometry used in analysis of lunar hydromagnetic shock wave 

and post shock atmosphere. 


The jump equations (Eqs. (17) to (21)) can now be simplified to 

[jOuJ = 0 

[u y ] =0 

[u z ] = 0 


and 



+ P + 



= 0 



(23) 

(24) 

(25) 

(26) 
(27) 


Taking the curl of the Ohms Law (Eq. 9) and solving for the jump equation yields 

[B z u x ]=0 < 28) 

o 
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Also the energy equation yields the following jump equation 


^ 1 


kT 


u 


2 W + 


P 


B 2 

Z 

PP-n 


- 0 


( 29 ) 


where it is assumed that VP x VB = 0. Assuming that 2 W » 1/2 n kT , 

6 Z ^ 3 &./p 

i.e., the ratio of specific heats is two, and that the solar wind Mach number is 
very large compared to unity, the following limiting solutions are obtained: 


B z2 = 3 B 


z 1 


P 2 - 3 p x 

and 

u - u , (30 a,b,c,d) 

y2 yl 

The hydromagnetic shock wave standoff distance A , is approximately given by 
the ordinary fluid shock standoff distance which is 


R 0 (*P\/ p 2) 

1 + /2 P l /~P 2 


(31) 


where R 0 is the radius of the moon. The maximum value of the stagnation point 
density calculated assuming P cc p 2 is approximately given by 


(32) 

This result indicates that radial variations near the stagnation point are quite 
small. However, it should be pointed out that, both the standoff distance and 
stagnation point density, as calculated, use the assumption that the net gas flux 
at the surface of the sphere is zero. If the sphere is losing mass, the net surface 
flux is not zero and both the surface density and standoff distance are increased. 
The ratio of the mass flux through the shock to the net surface flux is orders of 
magnitude less than unity in the case of the moon (Hinton and Taeusch (1964)) 
however, for example, in the case of comets, this ratio can be much greater than 
one. 
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The Post-Shock Lunar Atmosphere 


The individual particle densities and flux densities behind the shock wave 
can now be determined. Because of the formation of the shock wave, no slow 
neutral atoms or slow hydrogen ions can be present upstream of the shock, 
otherwise the free stream flow would be perturbed ahead of the shock by ioniza- 
tion of the neutral atoms and a stationary shock would not exist. Therefore the 
pre-shock flux is composed of only solar wind protons and electrons. Then the 
particle density behind the shock is 


3n *i 


n t 2 + n i2 + <2 + n a2 


(33) 


assuming that n e ~ n. , m i = m a and neglecting terms the order of m e /m r Since 
no charge exchange occurs ahead of the shock or in the shock, the energetic 
neutral density behind the shock is zero, i.e., n* 2 = 0. Then (33) in terms of 
the total ion density is 


3n u =n i2 + n a2 


(34) 


The mass flux density through the shock is conserved, therefore 


A - n * v* - n* V* + n s V s + n s V s 

V 1 “ "il 'll i2 V 12 + 12 12 + a2 V a2 


(35) 


where, again, terms the order of m e /m i are neglected V is the particle velocity 
normal to the shock and 4> i is the solar wind flux density. 

At the surface of the moon the total mass flux is zero or 


^o= n aoU°o = - n ioU i 0 -n: o u: o 


(36) 


where <K 0 is the slow ion flux density and U is the radial velocity in the plane 
4> =tt /2 (See Figure 2). The velocity U® 0 is the thermal velocity of the hydro- 
gen atoms emitted from the lunar surface. Since (1) U a0 , for a temperature of 
400°K, is much smaller than the solar wind velocity and, (2) the slow ion density 
at the surface is zero, the following relation is obtained: 



* 


iO 


Uto KoKo 


» 1 


u s 


aO 


n i0 U a0 


(37) 
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The total particle density at the lunar surface n Q is then given by 


which is the boundary condition on n* 0 at the stagnation point since n 0 is known 
from Equation (32). Therefore, the stagnation point slow neutral particle flux 
density is orders of magnitude less than the solar wind flux and, from Equation 
(35), the total ion flux density behind the shock is approximately equal to the solar 
wind flux density or 


n i2 V i2 * *1 


(39) 


The n. 2 and V i2 are 


and 


n i2 


= 3 n* — n 
1 1 


S 

a2 


V. = 1 (40 a, b) 

3 n . . - n s „ 

ll a 2 

where the only unknown is n® 2 . However, since n® 0 is known, it is possible to 
obtain~n® 2 from a solution of the equation of motion and the continuity equation 
for the slow ions . 

A detailed solution for the motion of slow ions behind the shock is beyond 
the scope of this paper, however, an approximate solution can easily be obtained 
if it is assumed that the slow neutrals are emitted from the surface at the thermal 
velocity in only the radial direction. Further, since collisions do not occur, the 
radial velocity is assumed to be relatively constant, i.e., U* ~ U® 0 . The con- 
tinuity equation for the slow neutrals then becomes 


1 JL 

r 2 Br 


(r 2 n®) = - 


n 


U s 

a 0 


( v i 


<t>* cr ’ 
1 cx J 


(41) 


where $* is the energetic proton flux density which is assumed to be constant, 
cr cx is the charge exchange cross section and v. is the photoionization rate. 
The solution to (41) is 


n r 


n 


(r 2 /R 2 0 ) 


exp. 


Or 
l cx 


u s 

U a0 


) R 0 


(1 - r/R 0 )> 


(42) 
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and the ion particle density for <f> = tt/2 and 8 - 0 is given by 


n. 

1 


S 

A; Il/v fl 

0 a 


( 43 ) 


if it is assumed that p H p 2 behind the shock. The neutral hydrogen and proton 
densities as a function of radius are shown in Figure 3 for v* = 5 x 10 _7 /sec, 
cr cx = 2 x 10" 15 cm 2 , =? 10 9 /cm 2 sec, n 0 = 3 x 10/cm 3 and Ug 0 = 4 x 10 s 

cm/sec (Hinton and Taeusch (1964)). The 9 variation of n®, to a first approxi- 
mation goes as the cos 9. The normal ion velocity and n® 2 along the shock 
(Equation (40b)) as functions of 8 are shown in Figure 4. 


Electric Fields in the Lunar Atmosphere 


An estimate of the electric field in the vicinity of the stagnation point can 
be obtained by evaluating the generalized Ohms Law at 8 = 0 assuming pH p 2 
behind the shock. This yields 

E = - u x B , — L VP a (44) 

n e e 



Figure 3-The thermal hydrogen atom and ion number densities behind the shock wave as 
functions of radius for <$> ~tt/2 and (9 = 0 and n^ = 10/cm 3 . 
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Figure 4-Ion velocity and thermal hydrogen atom density immediately behind the shock 
wave as a function of 9 , for 0 = tt/2 n u = 10/cm 3 0 X = 10 9 /cm 2 sec. Shock wave position 
is R g ^ (R q + A)/cos 6 . 


The electron partial pressure p e , from Equation (14), can be expressed as 



assuming that B/n e and W jVB are constants (Chandrasekhar (I960)). Then, since 
n e ^ n i and 



form Equation (27), the electric field is 
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(46) 


E = - u x B — 


B 2 >0 P el 

B? 


- 4 


Vn. 


M 0 e ( n i2 ) 2 

The mass center fluid velocity u, for 4> = tt/2 and 0 equal to zero, has only a 
radial component which is. approximately given by 


u r ~ 3 u x2 (1 - Rg/r 2 ) 


(47) 


The radial and tangential values of the electric field are 


%S3u xlBl (1 -R 2 /r 2 ) 


and 


2n 0 B 2 (O&L-A 


B i 


^ R o V 


M 0 e (n. 2 ) 2 R 0 


(48 a,b) 


at 6 = 0 in the <£ = tt/2 plane. The orders of magnitude of E^ and E r are 

E^ - 0 (u xl Bj) ~ 10~ 2 volts/meter 


and 


E r = 0 


"o B ? 


^0 e ( n i2^ 2 R 0 


£ 10“ 5 to 10 -4 volts/meter. 


Plots of E r and E^ are presented in Figure 5. 

The influence of the lunar gravitational field on the motion of the ionized 
lunar atmosphere and on the atmospheric electric fields is extremely slight. 
Inclusion of gravitational forces adds the following term to Equation (44) 



which has a maximum possible value of 10 -11 volts/meter at the lunar surface. 
Elsewhere in the atmosphere the space charge is quite small, therefore, E grav 
is essentially zero. 


58 




Figure 5-The radial and <p components of the lunar atmospheric electric fields as 
a function of radius for <P = tt/2 and d - 0. 

At the lunar surface, photo emission and secondary emission gives rise to 
an outward flux of electrons , most of which must be reflected back to the surface 
by an electric field. This radial electric field is provided by the net positive 
surface charge due to the decrease of surface electrons which are involved in 
photo emission. Next to the surface is an electron sheath which, as calculated 
by Singer and Walker (1962), is the order of several centimeters thick. To 
maintain the proper boundary conditions for the atmospheric electric field, which 
is radially inward, the sheath contains a surface density of electrons in excess 
of those needed to balance the positive surface charge. The electric fields 
within ten centimeters of the surface are the order of 10 volts /meter and de- 
crease very rapidly to the value of 10 ~ 4 volts/meter a few meters from the 
surface. 


Conclusions 

It is very probable that a shock wave of some kind is present on the moon 
and measurements by Ness (1966) have tentatively verified the existence of a 
lunar wake which might have been generated by a shock wave. The semi- 
cylindrical geometry of the shock, which has been proposed, results from the 
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condition that the magnetic field normal to shock must be zero for a steady-state 
shock wave. The steady-state semi -cylindrical shock wave cannot, of course, 
extend indefinitely along the magnetic field lines therefore, a non-steady shock 
structure must exist on each side of the steady shock. 

The application of the hydromagnetic jump conditions across the shock 
yields an increase in the total density and magnetic field of approximately three 
times the solar wind values. The stand-off distance of the shock is the order of 
0.2 times the lunar radius. The collisionless hydromagnetic model does not 
provide a description of the internal structure of the shock. However, it can 
be concluded that some form of enhanced ionization process, such as the impact 
ionization suggested by Beard (1966), operates in the shock because no neutrals 
can be present upstream of the shock. The only way in which neutrals arriving 
at the shock from the lunar surface can be reflected by the shock is by ioniza- 
tion and acceleration back into the flow behind the shock. The thickness of the 
shock is probably related to the ionization processes in the shock. 

The primary conclusions which are reached from this analysis are: (1) that 
most of the ions on the lunar atmosphere are swept along with the solar wind and 
interplanetary magnetic field, therefore, the lunar atmosphere is extremely tenuous, 
having a total particle density of no greater than 50/cm 3 , and (2) the electric 
field in the lunar atmosphere has a tangential component proportional to u x B 
which is the order of 10" 2 volts/meter and a radial component, proportional to 
the gradient of the electric partial pressure, which is the order of 10~ 4 volts/ 
meter. 
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A STUDY OF THE FREQUENCY AND STABILITY 
CHARACTERISTICS OF A C0 2 LASER 

by 

B. J. Graham* 


INTRODUCTION 

For deep space communication applications the C0 2 laser must be stable 
with respect to power and frequency output. A study of the output of existing 
lasers has been made to determine the design parameters of C0 2 lasers which 
will give a stable output. 


LASER DESIGNS 

Figure 1 is a schematic of the first laser studied. This laser had internal 
gold mirrors sealed by a bellows arrangement for external adjustment within 
a 1.2 meter tube of 1-3/16 inch inside diameter. The mirrors at each end of the 
tube were nearly 100% reflecting with the exception of a small hole (various 
mirrors had holes with diameters from 1 mm to 5 mm) in the output mirror. 

The system was basically a flow system of C0 2 , N 2 , He and many combinations 
of partial pressures of these gases were tried. A glass jacket around the dis- 
charge tube was used for air cooling. The laser, rf excited by a Johnson Viking 
80 watt transmitter, had an output of approximately 3 watts. A second laser 
very similar to the first was dc excited by electrodes in each end. The external 
jacket in this case was used for either air or water cooling. 

From the spectral data taken from these lasers a second basic design 
evolved (Figure 2). 

This laser had a discharge tube sealed at both ends by salt windows set at 
Brewster's angles, and with the mirrors externally mounted. The tube was only 
30 cm long and contained an external glass jacket for cooling by circulation of 
air or a liquid dielectric. The mirrors were readily interchangeable and various 
configurations of radius of curvature of the mirrors and size of the output holes 
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Electrodes for 






Loose bushing to 
allow free expansion 



were tried. All of these lasers were rf excited but the power output was low (50 
milliwatts) compared to the longer lasers. Rf excitation was accomplished by 
placing parallel electrodes along the sides of the laser tube. The electrodes 
and plasma tube formed a capacitive system which was tuned to the resonant 
frequency of the rf exciter by addition of an inductance coil and adjustable tuning 
capacitor. 


SPECTRAL CHARACTERISTICS 

The laser action for all the lasers used was from either the P- or R- branch 
rotational transitions of the 00°1 - 10°0 vibrational bands of C0 2 (Reference 1). 
The spectral output of the lasers either rf or dc excited, when the mirrors were 
sealed within the tube, was essentially the same and was quite unstable. All 
observations where made with the Cary White 90 infrared spectrophotometer. 

The laser action varied from a single transition to as many as six transi- 
tions simultaneously but these were seldom repeatable. Figures 3, 4, and 5 are 
typical of many spectra taken of both the rf and dc excited lasers with sealed 
windows. On some runs spectral lines were separated by as much as 36 cm -1 . 

It was determined that the fluctuations of input power, tube pressure, and thermal 
effects all contributed to the instabilities. 

The shorter lasers with external, isolated mirrors did not support as many 
different spectral lines. In almost every run frequency stability had greatly 
improved over the previous lasers . As many as four transitions were observed 
at once, but all transitions were tightly grouped with the spread between lines 
never greater than 2 cm - 1 . Figure 6 shows a very sharp single spectral line 
from one of the short lasers. Figure 7 is the same spectral line but under a 
higher resolution and dispersion. The peaks shown here are much too close to 
be spectral lines. They are probably some of the more dominant modes within 
the P(20) line brought out by the higher resolution of the spectrophotometer. 


MODES 

The laser cavity supports not only the spectral radiation, but the electro- 
magnetic modes which are characteristic of the resonant cavity. For a confocal 
resonator the resonant condition is given by (Reference 2): 

4b/X. = 2q + (m +n + 1) 
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P ( 16) - 10,55 y 



Figure 4— Spectrum of the Radiation From a D.C. Excited C0 2 La 
Containing C0 2 Only 
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Figure 6- Very Sharp Laser Line (P(20)) Under 
Moderate Resolution. This Line Can, With Some 
Adjustment, Be Regularly Obtained From The 
Short C0 2 Lasers With the Isolated Mirrors. 
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where b is the mirror separation, X the wavelength, q the number of half-waves 
along the length of the laser, and m and n designate the mode. For a laser 
resonator of 50 cm length the separation of the electromagnetic modes is of the 
order of 5 x 10 ~ 3 cm -1 . This means that as many as 20 electromagnetic modes 
can appear under the width of a single spectral line (natural line width ~0.10 
cm -1 ). An experiment was designed (Figure 8) to look at these individual modes 
within a spectral line. In the experiment, a Lansing scanning interferometer 
was modified by replacing the flat Fabry-Perot type mirrors with a confocal 
arrangement of nera 100% reflecting gold mirrors, each with a 0.25 mm diameter 
hole drilled in its center. Radiation is coupled into the cavity through the hole 
in one mirror. When the cavity mirrors are critically aligned, for a given 
wavelength, i.e. the wavelength of an electromagnetic mode, the cavity will be 
in resonance and radiation from the output hole will be increased. When the 
length of the cavity is changed, this wavelength will not resonate and the output 
energy is reduced. Part of the incident radiation is picked off by a beam splitter 
and sent into an infrared spectrophotometer. This will make it possible to de- 
termine the wavelength of the spectral line from the laser, although it cannot 
clearly resolve the modes found within a spectral line. 

One of the mirrors of the Lansing interferometer is backed by a piezoelectric 
translator which is driven by a linear ramp generator. The total translation in 
this particular arrangement is 1.2 microns which is a little more than one-tenth 
of the incident wavelength (10.6 microns) from the CO a laser. 

Figure 9a shows an oscilloscope sweep of the laser signal coupled out of the 
interferometer when the mirrors are stationary. Only the instrument noise 
appears in this trace. However, Figure 9b shows a sweep of the signal as the 
length of the cavity is shortened by 1.2 microns. The spikes on this trace are 
the electromagnetic modes that are brought into resonance as the cavity length 
is changed. , 

At this writing there has not been a detailed analysis of these traces. This 
could be a tedious task because the spikes along the trace not only represent the 
modes reinforced within the cavity but also the degenerate modes. 

Figure 10a shows a sweep of the modes, and Figure 10b shows a sweep after 
the laser mirrors have been adjusted very slightly. This figure shows additional 
modes being reinforced and strengthened and could thus give an indication of 
what goes on inside the laser itself. 

This interferometer arrangement should make it possible to completely 
analyze the closely spaced electromagnetic modes of the C0 2 laser. It is hoped 
that this knowledge will aid in the design and development of a highly stable laser. 
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Figure 9(a)-Signal Plus Noise From CO 2 Laser Radiation Through 
the Apparatus in Figure 8. (b)-Spikes in the Sweep Represent 
Laser Mo'des as They Brought Into Resonance as the Resonator 
Cavity is Shortened. 
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Figure 10(a)— Spike.s in the Sweep are Laser Modes as They are 
Brought Into Resonance By Shortening the Resonator Cavity (Fig- 
ure 8). (b)— The End Mirrors of the Laser are Slightly Adjusted 
and Some Modes are then Reinforced Over Others. 
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CONCLUSIONS 


The design efforts on the short lasers with external, isolated mirrors were 
successful in producing a more monochromatic source. Study of this design 
shows that more needs to be understood about the following areas: 

a. excitation mechanism of the laser action, 

b. the modes within the laser cavity and how to control them, 

c. the dependence of the lifetime of a closed laser tube upon the gases used 
and the preparation of the laser tube. 

The experimental set up in Figure 8 is the first time a cavity resonator has 
been used in the infrared wavelength region. This has much promise as a means 
of analysis of the modes of lasers in the infrared region. This technique has good 
possibilities of application in the visible radiation region also. 


REFERENCES 

1. C. K. N. Patel, "Continuous — Wave Laser Action on Vibrational — Rotational 
Transitions of C0 2 ," Physical Review, Vol. 136, No. 5A, 30 Nov. 1964, pp. 
A1187 - A1193. 

2. G. D Boyd and J. P. Gordon, "Confocal Multimode Resonator for Millimeter 
Through Optical Wavelength Masers," The Bell System Technical Journal, 
March 1961, pp. 489-505. 

3. T. K. McCubbin, Jr., Ronald Darone, and James Sorrell, "Determination of 
Vibration — Rotation Line Strengths and Widths Using a C0 2 - N 2 Laser," 
Applied Physics Letters, to be published. 





RELATIVE PHASE DIFFERENCE ANALYSIS PROGRAM 

by 

D. M. Koenig* 


INTRODUCTION 

Oscillator-pair data, which is relative phase difference ( 86 ), is recorded 
continuously on strip charts (see Figure 1). 



Figure 1 


This curve is sampled at times t. , i ■= 1, 2, . . N 3 where L = (i - l)At 
with At usually one hour. The temperature at these times is also sampled; thus 
a set of N data triples is produced: (T. , S0. , t. ), i = 1, 2, . . .N . 

The data are discontinuous because of scale shifts and crossovers. These 
points of discontinuity are coded (values are assigned to the variable L . ) ap- 
propriately so each sample point now has four numbers associated it: t . , T. , 
S6 l i , L i . For every data point, T , S<9. , L £ are punched onto one IBM card. 

These cards are read by the program and the data is corrected for the 
discontinuities. The corrected data is printed out for checking purposes. The 
second phase consists of fitting the data to a double polynomial in time and 
temperature, i.e., the coefficients A k are found in 

'*rN t V. t V! , A, T i t ¥! ,e i 
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by minimizing 


N 



i = 1 


This is done by a standard least squares procedure which requires the solution 
of a 5X5 linear system. 

The estimated standard deviation s , where 


N 



i = 1 


is calculated and each e. is compared to 3s . If |e. | > 3s , the ith point is 
rejected. After all of the e. are compared to 3s the coefficients A k are re- 
calculated using only those points passing the rejection test. 

The last phase consists of an autocorrelation of the residuals e. to show 
the presence of any periodic behavior. 


DETAILED DESCRIPTION OF SUBPROGRAMS 
A. Overall Configuration 

The program consists of seven subprograms entitled as PRG1, PRG2, etc. 
The overall flow chart is shown in Figure 2. The statement ’’transfer to PRG1” 
means transfer to the appropriate routine, carry out the necessary operations 
and then return to the point of transfer. 


B. PRG1: Data Corrector Routine 
1. Critical Point Coding 

Considering the case of increasing data, discontinuities of three types can 
occur and these are illustrated below in Figure 4. 
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PRG3: COUNT GOOD 
AND BAD PTS. 



Figure 2 
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<G 2: | CALCULATE AVERAGE 
TEMPERATURE AND SUBTRACT 
FROM ALL TEMPERATURES 


RETURN 


RG 4 : j REDUCE THE MATRIX BY 
THE GANSS-JORDAN METHOD AND 
DETERMINE THE COEFFICIENTS A|< 


RETURN 


PRG 5: 1 CALCULATE THE e*, , e , e 2 
AND SUBTRACT e FROM EACH ej 


RETURN 


JLj CALCULATE s, COMPARE ej 
TO 3s AND TEST FOR REJECTION 


RETURN 


PRG 7: 


AUTO CORRELATE THE e ; AND 


PRINT OUT THE AUTO CORRELATION 


Figure 3 
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CROSSOVER UPSHIFT 


DOWNSHIFT 


PHASE OR 
FREQUENCY 



Figure 4 


Often one or more points are missing during a shift. To account for this, 
a certain value is given to the code variable which is called LABJ, where J re- 
fers to the Jth point. The coding system is given in Table 1. 


Table 1 


Type of Point 

Code 

last good point before shift 

1 

bad or missing point 

3 

first good point after crossover 

5 

the very last good point 

6 


An example of the point coding is shown in Figure 5. If a point does not fall 
into any of the classifications of Table 1 it is left uncoded. Thus for each data 
point one IBM card is prepared containing a temperature, relative phase dif- 
ference, and a code (if it is a critical point). 


2. Overall Procedure for PRG1 

PRG1 reads and processes the data cards one at a time. On reading a card, 
and thus a data triple, values are temporarily assigned to the variables TMP, 
RPD, LABJ which stand for temperature, relative phase difference, and label. 
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5 1 3 


6 LABJ 



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 J 

Figure 5 


Depending on LABJ, a new value is assigned to RPDB which is the corrected 
relative phase difference. The values of TMP, RPDB, and LABJ are assigned 
to DTI(J), OMEGA(J), and LAB(J) where J is the point index. Having stored the 
corrected data, another card is read in , new values are assigned to TMP, 
RPD, LABJ, J is stepped up by one, and the cycle is repeated. 


3. Correction for Crossovers 

Each point is corrected by adding to it the value of a variable called AOF 
or "add-on-factor." Initially AOF is given the value zero. After each crossover 
AOF is increased by 1000. Since the data triples (TMP, RPD, LABJ) are read 
in sequentially, AOF is changed whenever LABJ has the value 5. 


4. Correction for Shifts 

If a point's code is 1, it is the last good point before a shift. To correct 
the next good point a straight line is passed through the previous data points 
back to a shift or a crossover. This straight line is used to estimate the value 
of relative phase difference at the next good point had there been no shift. The 
difference between the actual and the estimated relative phase difference is 
subtracted from AOF and is used to correct the current point. The value of 
AOF is left unmodified until another crossover or shift occurs. To calculate 
the intercept and slope of the straight line a least squares method is used. 
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5. Modification for Decreasing Data 

For decreasing data the add-on-factor becomes a subtract-off -factor. To 
account for this, AOF is multiplied by a factor FAC which is either plus or minus 
one. The value of FAC is read in at the beginning of the program and is set to 
plus one for increasing data. 


6. Flow Diagram for PRG1 
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C. PRG2 : Subtraction of Mean Temperature 

The mean temperature is calculated and then subtracted from each individual 
temperature. 

D. PRG3: Polynomial Fitting Routine 

The primary task of this program is to determine the coefficients A k in the 
double polynomial. The method of least squares is used and is derived as 
follows : 


se { = Aj + A 2 T i + A,T* + V, * A s tf + e, 


=o 

^4- 1 


v -1 Be. 

) e. — i =0 

^ 1 


k = 1, 2, 3, 4, 5 


e. = SS. - A, - A 2Ti - - A 4 t. - ^ 


Be. 

=- !. 

BA 1 


' Be. 


BA, 


1 - - T. , 


Be. 

L - _ t2 

ba 3 


Be. 


BA 4 ** 


— t = - 1? 
ba 5 


2^89. = NAj + + A $ 22 t? + A + AsJ]'" 


22s9 1 T i = A^T, + A ,Jjl + A 3 22t? + A^tJ, + A^tfT, 
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These five equations are solved for the coefficients by means of a Gauss Jordan 
reduction routine (PRG4). 


E. PRG4: Gauss -Jordan Reduction Routine 

The procedure by which the linear system, constructed in PRG3, is solved 
is standard (see p. 224), Southworth and DeLeeuw, "Digital Computation and 
Numerical Methods," McGraw-Hill, 1965), and will not be discussed here. The 
title of the subroutine is GSRDN (A, B, X, N, DET). A is a doubly subscripted 
variable representing the elements of the coefficient matrix, B is a singly sub- 
scripted variable representing the input matrix and X is another singly sub- 
scripted variable representing the unknown matrix, 

(A) {X} = {B} . 

N is the dimension of (A) which for the case discussed in D is 5. DET is the 
value of the determinant of (A). 

F. PRG5: Calculation of Error at Each Point 

At each point e £ , i = 1, . . . n is calculated, where 

e. = 86. - Aj - A 2 T. - A 3 T2 - A 4 t. - AgtJ 

at those points for which LAB(I) is equal to 3, e. is set equal to zero. The 
average of the errors 
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i = l 


and the average of the errors squared 


n 



i= 1 


is calculated. 


G. PRG6: Rejection of Bad Points 

The standard deviation of the e. is calculated 




-i /&-(&)■/■ 

and then all the |e. | are compared to 3s . If | e i j is greater than 3s the value 
of e. is set to zero and LAB(I) to 3. The procedure is repeated using the re- 
maining non-zero e i where n is the number of non-zero The flow diagram 
for this routine is shown below. 
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AN INVESTIGATION OF METHODS FOR THE EXPERIMENTAL 
DETERMINATION OF PCM SIGNAL TO NOISE RATIO 

by 

Jack F. Morris 


INTRODUCTION 

The percentage of recovery of digital PCM data from a given magnetic tape 
can in general be related to the signal-to-noise ratio (SNR) at the telemetry 
system receiver. Occasionally, however, data that might normally be recovered 
is lost through processing errors or equipment malfunction. A study was pro- 
posed by the Data Processing Branch, the object of which was twofold: the 
telemetry link itself was to be investigated for the sources and nature of the noise 
at each stage of the data path through the receiver system to the analog data tape, 
and secondly, since the probability of bit error is directly related to the analog 
data SNR, it was proposed that one or more methods might be developed by which 
an estimate of data recovery might be made, independent of the processing line. 

Spectral characteristics of the receiver system noise and SNR's throughout 
the system were studied experimentally using the ground stations set up in the 
Information Processing Division (IPD) , and in the Network Test Training Facility 
(NTTF). The effect on data recovery of varying the modulation index and of 
varying the IF and post-detection bandwidths were also investigated. 

In the second phase of the work several of the classical methods of signal 
analysis were investigated for applicability to PCM analog data. Computer 
programs and digital techniques were developed for obtaining signal amplitude 
probability distributions over any desired segment of analog data, for obtaining 
bit-by-bit statistics within the frame sync pattern, and for obtaining a spectral 
analysis of the data. Each of these programs will yield a SNR for the analog 
data, based upon the statistics developed in the program. An estimate of the 
probability of bit error, and therefore an estimate of the quality of the data, or 
the recovery expected, can thus be made for any given PCM analog data rape. 


*Electrical Engineering Department University of Missouri at Roila. 



PHASE I - THE TELEMETRY SYSTEM 
SUMMARY OF THE WORK 

In an attempt to identify the sources of noise and the nature of the noise 
in the PCM telemetry receiver system, the following approaches were taken: 

1. Magnetic tape recorded data, in the OGO-D satellite format were trans- 
mitted through the collimating antenna to the Satan antenna at NTTF. 

The received data were recorded at a series of known signal levels, 
calculated from a reference receiver system noise figure. 

2. A clock signal, at a 64 kilobit persecond rate, was transmitted, received, 
and recorded in the same manner. The purpose of this was to provide a 
known signal, with a relatively simple frequency spectrum, to facilitate 
noise spectrum analysis. 

3. Using the receiver system in IPD, simulated data were inserted into 
the RF path of the system. The modulation index and the signal level 
were both fully controllable, and the effect on the output data due to 
various inputs was investigated. For low signal conditions an apparent 
decrease in SNR improvement through the PM demodulator was noted, 
as was expected, and for high signal levels the SNR improvement 
approached the theoretical value. 

4. In an investigation of the effect of tape recording on data quality, simu- 
lated data of known quality were recorded from the output of the IPD 
receiver system. The tapes were then played back into the IPD data 
analysis equipment. In each case few, if any, more errors were observed 
on playback than when recording, and the maximum increases in error 
rate represented possibly a 0.1 db degradation in SNR through recording. 


RESULTS OF RECEIVER SYSTEM INVESTIGATION 

The primary, and surprising, result obtained in working with the trans- 
mission of simulated data through the system at NTTF is that data quality is 
not always reliably represented by the signal level indication on the system 
receiver. Clean simulated data were recorded at 60 inches per second (ips) 
by IPD people. These tapes were at data rates of 1, 8, and 64 kilobits per 
second (kbps). The simulated data were transmitted through the collimating 
antenna, received by the system’s Satan antenna, and recorded on magnetic 
tape. The OGO-D satellite Operations Plan specifications were used for receiver 
and recorder settings. Four minutes of data were recorded at each of the signal 
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levels calculated for SNR as indicated in Table I below. The recorded data 
were then played back through the signal analysis line in the Processor Devel- 
opment Branch, IPD. A non-uniform discrepancy was noted to exist between 
calculated SNR (and, therefore, bit error rate) and the observed bit error rate, 
as indicated: 


Table I 

Calculated and Measured SNR 


SNR Calculated 

SNR measured 

Error 

by NTTF, in db 

by IPD, in db 

in db 

4.85 

3.00 

1.85 

7.30 

4.90 

2.40 

9.65 

6.80 

2.85 

11.90 

8.60 

3.30 

13.90 

10.20 

3.70 


Several possible explanations exist for the discrepancies noted, none of 
which has been fully investigated at this time. One of these possibilities is in 
the nature of the simulated data: because of the bandwidth limitations of mag- 
netic tape recording, it is difficult to obtain a satisfactory simulation of satellite 
data by recording and playing back the output of a signal generator. Another 
problem area is in the calculation of SNR using the IF bandwidth and system noise 
figure as the basis. Experimentally measured noise power in the IF seems to 
be slightly higher than the value calculated using the noise figure. It has been 
proposed that perhaps the IF filter cutoff is not sufficiently sharp that simply 
using the half-power bandwidth will determine total IF noise power to the desired 
accuracy. Further investigation of each of these problems is planned. 

A spectral analysis of the recorded clock signal plus noise revealed that 
the sky noise plus receiver noise in the system does have the "white" or uniform 
frequency distribution so often assumed in the literature pertaining to signal 
analysis. Observation on a spectrum analyzer, and measurement of noise power 
using various ranges of post-detection filter bandwidth, established the uniform 
nature of the noise in the 136 MHz range utilized by the telemetry system. 

Using the receiver system in IPD, experimental values of SNR improvement 
through the PM demodulator were obtained, along with results indicating essen- 
tially negligible degradation of the data through tape recording. The absolute 
noise power in the receiver IF measured in the absence of any signal and with 
manual control of the AGC. Signal levels were then adjusted for various SNR’s, 
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referred to the information bandwidth, and a frame sync pattern bit error count 
was made. Results were consistent over many runs and several days, and the 
values obtained are shown in Table II below: 


Table II 

SNR Improvement through PM Demodulator 


SNR in 

Post-detection 

SNR 

IF Filter 

SNR 

Improvement 

1.8 db 

2.5 db 

0.7 db 

4.8 

6.1 

1.3 

7.8 

9.6 

1.8 

10.8 

12.7 

1.9 


A modulation index of 1.25 was used for this work, which would yield a 
SNR improvement of 1.9 db for large signals, according to the theory. Although 
the last two measurements may be suspect, because of the very large number 
of frames of data which must be processed for valid statistical results, it should 
be noted that throughout all of this investigation the results obtained for data with 
SNR greater than about 9 db agreed in every case with theoretical predictions. 

To complete the receiver system phase of this initial investigation, the 
signals described above were recorded on magnetic tape and then played back 
to determine the amount of degradation of SNR in the recording process. For 
each recording the IF was six times the information rate and post detection 
bandwidth was ten times the information rate. Tape speed was 30 ips. In the 
worst case the degradation amounts to about 0.1 db, which is essentially negli- 
gible. Results appear in Table III. 


Table III 

Comparison of Error Rates 


SNB 

(db) 

Pre-recording 
error rate 
(bit errors/frames) 

Post-recording 
error rate 
(bit errors/frames) 

2.4 

1300/27,000 

1350/27,000 

2.5 

1340/27,000 

1350/27,000 

6.1 

141/27,000 

163/27,000 

6.1 

165/27,000 

165/27,000 

9.6 

14/100,000 

14/100,000 

9.6 

12/100,000 

14/100,000 
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Several more data points were taken, at different post detection bandwidths, 
and the results were similar to the above. A slight improvement was observed 
in some cases and a degradation in other cases when operating the system with 
narrowed IF and post detection bandwidths. It appears at this writing that a 
slight improvement in data recovery can be obtained, using narrow band filter 
settings. The improvement is very small, however, and is very much a function 
of the tuning and adjustment of the receiver. When the receiver was timed off 
of the carrier frequency by as little as 10 kHz the degradation in performance 
began to be noticeable: 200 errors in 27,000 frames, as opposed to 165 errors 
for the wider IF bandwidth data in the table above. The trade off between the 
system accuracy required to obtain a very slight improvement and the minor 
degradation experienced using the wider band system settings definitely favors 
the latter. 


PHASE II - PCM ANALOG DATA ANALYSIS 
STATEMENT OF THE PROBLEM 

Satellite telemetry data recorded on instrumentation type analog magnetic 
tape at the worldwide network of STADAN stations are mailed to Goddard Space 
Flight Center (GSFC) for processing. The telemetry data acquired from the 
satellites by the ground stations are received in less than an ideal form that 
includes blankouts, weak signals and noise. Phase II of this investigation was 
concerned with the study of PCM signal and noise characteristics of actual 
flight data recorded on magnetic tape, and the development of computer tech- 
niques and procedures for automatic signal evaluation. 

Most of the theoretical investigations of PCM SNR which have appeared in 
the literature have been based upon an ideal square wave type signal corrupted 
by white Gaussian amplitude noise. Under these assumptions, it is easily shown 
that the square of the mean value of the signal-plus -noise represents the signal 
power and the variance about the mean represents the noise power. That is: 

mean squared m 2 signal power 

^ — = = SNR 

variance a ' 1 noise power 

Thus, when experimentally measured values of signal mean and variance 
are used as an estimator of SNR, in the statistical sense, there are standard 
statistical procedures by which the true SNR may be obtained with any level of 
confidence. All that is required by these procedures is that a larger and larger 
number of samples be used to obtain the m 2 / o- 2 estimator of NSR as the confi- 
dence level is raised on the true value. 
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When the actual flight data signal-plus-noise is considered, however, it is 
seen that the signal is not ideal and the noise is not strictly an amplitude noise. 
Rather, the signal portion of the waveform represents an ideal square wave 
which may have been pre-filtered before transmission, has been subjected to. 
atmospheric fading and phase shift, has been further filtered by the IF and post- 
detection filters of the receiver system, and again filtered by the transfer char- 
acteristics of the record qnd playback electronics of the tape system. A typical 
flight data waveform is represented by Figure 1 below: 



Figure 1-Typical PCM flight data waveform. 


The filtered signal waveform is in addition corrupted by amplitude noise, 
by phase noise or phase jitter, by wow and flutter of the tape recording system 
and by other influences which may be considered to be noise, such as the spin 
modulation introduced by the rotation of the satellite. 

The processing of PCM analog data is further complicated by the fact that 
it must generate its own timing signals, or clock. The data tape is played into 
a bit synchronizer and signal conditioner which must lock onto the data bit stream 
with an internal voltage controlled oscillator (VCO). The VCO frequency is 
corrected automatically for small changes in signal frequency, but cannot readily 
synchronize with data of low SNR. Thus, for the low SNR condition a reliable 
data-rate clock, derived from the VCO, may not be available. 
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SUMMARY OF THE WORK 


In an attempt to obtain a measurement of PCM data SNR, auto-correlation, 
cross-correlation, and ensemble averaging of the analog bit stream were 
considered. Except for an adaptation of the cross -correlation technique, from 
which a synchronizing signal was obtained, the corelation process did not lend 
itself to application to this problem. The reason for this lies in the nature of 
the data format, as explained in the following section. 

Synchronous and random sampling procedures were used to acquire statis- 
tics on the data from which signal amplitude probability distributions, signal 
mean and variance, and spectral characteristics of the signal were determined 
through digital techniques. The amplitude distribution program was designed 
to yield a SNR for the data from the mean squared signal, representing signal- 
plus-noise power, and the mode, representing the signal power alone. The signal 
mean and variance program was designed to yield average and/or bit-by-bit 
statistics, from which SNR was calculated directly. Agreement between the two 
independent programs was found to be consistent, although a calibration routine 
is required before agreement is complete. 

A spectral analysis program was written utilizing the same data samples 
as above and an adaptation of the Cooley-Tukey fast Fourier transform (FFT) 
algorithm. This program was designed to yield a SNR based upon the knowledge 
that the second harmonic of the signal bit rate is essentially non-existent. Any 
frequency components in the neighborhood of the second harmonic may therefore 
be interpreted as noise, and a noise power figure obtained. 


RESULTS OF PCM ANALOG SIGNAL INVESTIGATION 

The data format of the OGO-D satellite was the format used for all of the 
signal investigation. This format is typical of many PCM systems, and consists 
of a "frame" of 128 nine-bit "words" . Each frame has three known or programmed 
words used to locate the frame and synchronize the words for computer manipula- 
tion. These three words are termed the "frame sync pattern", and in the OGO 
format are the first three words of the frame, one shown schematically in 
Figure 2. 

Now consider that the analog signal may be represented by 

f (t) =S(t) + N ( t) , 


95 



FRAME 

SYNC 


123 NINE-BIT WORDS 


FRAME 

SYNC 


-v 


27 BITS 


1125 BITS 


Figure 2— PCM data format for OGO satellite, 

where S(t) is a periodical signal with zero mean, and N(t) is a random noise 
with zero mean and without any periodic component. Then the auto-correlation 
function of f(t) is given by 


0 ff (T) = [s(t) + N(t)] [S(t + T> + N(t + T )] 

- ^ss( T ) + <W T ) + <^s( r ) + ^sn( t ) 

“ 0ss( r ) + c hiN^ T ^ 

under the conditions of zero mean signal and noise. 

In a practical situation 4> m (r) tends to zero at some sufficiently large value 
of such that for uncorrelated signal and random noise the function 


4>fi (r) = <p ss (r) , for larger, 

as indicated schematically in Figure 3. From the defining equation it is seen 
that 4> {f (r - 0) represents the mean squared signal-plus -noise, or the total 
signal-plus-noise power. 0 ss (r = 0) represents the signal power alone, and its 
value may be determined at some t-t 0 such th.2A,4> m ('r 0 ) = 0 as indicated in 
Figure 3. 

One way in which the auto-correlation process might be applied to PCM 
data, in an attempt to determine the signal component and the noise component 
of the data waveform, would be to treat the known frame sync pattern as a 
periodic function in the presence of noise. To do this requires that the data be 
considered purely random, as well as the true noise component. The limited 
usefulness of the method becomes apparent when it is noted that the frame sync 
pattern represents only 3 of 128 words in the frame, or less than 3% of the total 
signal. Thus, even though measurement of 0 ss (r o ) may be practical, the errors 
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introduced by any periodic data would in all likelihood be greater than the 
measurement. That is, it would only require that 3 data words be constant for 
several consecutive frames to introduce a 100% error in the signal power 
measurement. 



Figure 3— Autocorrelation function of a periodic signal plus random noise. 


Another approach is to consider the frame sync pattern as the periodic 
component and the noise in the frame sync words as the random component. 
This requires gating out the data, and considering only the first 27 bits of each 
frame, as shown below: 


FRAME 


FRAME 


FRAME 

SYNC 

kk 1 

SYNC 

kk 

SYNC 


17ms 


1125 BITS 


TPbhs 


1125 BITS 


GATED OUT 


GATED OUT 


Figure 4— Schematic of conditions for auto-correlating the frame sync words. 


This latter technique was not pursued in the course of the investigation. 

The timing problem associated with gating out the frame sync patters and 
essentially connecting them end to end coherently was not solved during the 
time allotted. Further work might be done in this area of a manageable solution 
to the timing difficulty can be found. 
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The most useful technique discovered during this work was an adaptation 
of an amplitude probability distribution program. A procedure was developed 
under which either random samples or synchronous samples of the bit stream 
were collected, the probability of each amplitude increment calculated, and the 
probability as a function of amplitude plotted by digital computer. Figure 5 shows 
a representative probability curve for a data waveform such as that in Figure 1: 



The mean squared value of the distribution is proportional to the total 
signal-plus-noise power when the mean value is zero. That is, 


v 2 (t) = E (v 2 ) 



v 2 p(v) dv 


or E(v 2 ) = j^E(v) J + cr 2 , 

where a 2 is the variance, or second moment about the mean. When 

|E(v) ] 2 = 0, Then E(v 2 ) = a 2 

is the power delivered to a one ohm resistor across which v is the voltage. 

Under the assumption that the mode E of the distribution is related to the 
signal power, the SNR was expressed as 

signal power 

SNR = — — ^ — T r 

(signal + noise power) - (signal power) 
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or 


SNR = 


kE 2 


cr 2 - kE 2 , 


where k is a proportionality constant, determined by calibration. 

This program was used with samples of actual flight data for which error 
statistics had previously been obtained. Accurate agreement over a range of 
SNR's from approximately 4 db to 10 db was obtained using a value of k =0.78, 
where the constant was chosen to yield a SNR of 20 db for clean simulated data 
played back from magnetic tape. The program generated consistent results when 
used with random samples of data, synchronous samples of data, or signal 
generator waveforms, and is considered a successful technique for obtaining 
SNR’s. 

An oscilloscope display of the PCM flight data reveals a slight pattern 
sensitivity in the data waveform. That is, certain bit patterns in the data seem 
to introduce anomalies in the waveform that are not observed with different bit 
combinations. A Polaroid photograph of several consecutive sweeps through 
the same data word is reproduced in Figure 6. Several repeating high and low 
values of the waveform are clearly seen in this picture. 



Figure 6— Polaroid photo of PCM data waveform. 


To investigate the apparent pattern sensitivity of the waveform, a digital 
computer program was developed which uses a cross-correlation technique to 
locate the frame sync pattern in a stream of data samples obtained synchro- 
nously, one or more samples per information bit. When the frame sync pattern 
is found, the samples from each bit within tne pattern are stored separately. 
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The mean and variance about the mean for each bit in the frame sync can be 
computed, as well as the mean of the 27 bits in the pattern taken as a whole. 

Other statistics using these samples are easily generated with minor program 
changes. 

It was found from this program that the pattern sensitivity was clearly 
revealed by the samples. Figure 7 represents the mean value, by bit position, 
obtained for the 27 bits of the frame sync pattern on a typical flight data tape 
of approximately 12 db SNR. The 3 a value, where a is the standard deviation, 
is shown on each bit and represents the noise level that would results in an error 
rate of approximately 5 x 10~ 3 

In this program the SNR is derived from the ratio of mean squared to 
variance, where in this case the mean is determined at a constant phase point 
on each bit by synchronous sampling . The variance is taken about the signal 
bit mean, and is proportional to the noise power. It was found that the variance 
of the noise power from bit to bit is essentially negligible. That is, the noise 
contribution at each bit is the same, regardless of the mean value of the signal 
bit. This can be seen in the relatively uniform heights of the 3 o- points in Fig- 
ure 7. Thus for high noise conditions certain bits would tend to be in error more 
often than others. That this is indeed the case is borne out by bit error statis- 
tics developed independently of this program by IPD. 

A valuable feature of the amplitude distribution and cross -correlation 
programs discussed above is that they are applicable without modification to 
the study of the PCM bit synchronizer output SNR. All that is necessary is 
that the output of the ’’integrate and dump” circuitry be made available for use 
as a signal source. Amplitude distributions, mean and variance, and SNR 
figures for the output data may be compared over whole data tapes , or over any 
desired segment of data. This feature should contribute significantly to an 
evaluation of bit synchronizer performance. 

The last program developed during this work was an adaptation of the 
Cooley-Tukey fast Fouruer transform (FFT) algorithm. The program will yield 
a frequency spectrum with a maximum width corresponding to one half the 
sampling rate. The resolution is a function of the number of samples processed 
for any given spectral ’’display”. There is some danger of losing or distorting 
spectral information near the highest frequencies, and the sampling rate should 
therefore be chosen somewhat higher than twice the maximum frequency of 
interest. 
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When the PCM analog data is applied to a conventional swept frequency 
spectrum analyzer, or to a real time analyzer, it is seen that the second har- 
monic component of the signal is essentially zero. The spectral null near the 
second harmonic in quite broad, and since the signal waveform is ideally a square 
wave type signal, which has undergone linear filtering, it should be expected 
that this would be true. Taking advantage of the null in the signal spectrum, and 
assuming that the frequency components near the second harmonic are due to 
noise alone, a digital filter is placed around any desired band of frequencies 
near the signal null. The digital filter is used to find the noise power density, 
and therefore the total noise power in any desired bandwidth. The FFT program 
will also find the total signal plus noise power, and the SNR may be obtained for 
the block of data from which the samples were taken. 

Calibration of the digital noise filter is required for each different sample 
rate, data rate, or frequency resolution used in an analysis. The calibration is 
essentially a program input which determines the frequency width and location 
of the filter, and is easily determined from the data being analysed. Excellent 
agreement between results obtained with the programs previously discussed and 
the results of the FFT program has been observed. In addition, it should be 
noted that the FFT program has the capability of examining as small a segment 
of data as desired, which is of benefit in an evaluation of burst noise or other 
localized problems. For this reason the FFT is a useful addition to the library 
of data analysis programs. 


DISCUSSION OF FUTURE WORK 

Calibration and refinement of the digital techniques developed during this 
initial investigation will be continued, using simulated data and flight data with 
known bit error statistics. 

Additional planned future work includes further analytical and experimental 
study of the telemetry receiving system, the analog tape recorder, and the 
analog to digital processing line. A digital model of the bit synchronizer and 
signal conditioner will be developed and studied. An all digital model of the 
data gathering and reduction process might be considered the goal. 
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COMPUTER SYSTEMS FOR 
REAL-TIME APPLICATIONS 

by 

David F. Nyman* 


ABSTRACT 


The STADAN tracking sites are currently operating with very 
little automatic control of the equipment involved with tracking a 
satellite. The Advanced Development Division at Goddard Flight 
Center is currently studyng the possibility of using a computer to 
control the tracking station equipment. This report will describe 
the functions of a computer and its associated monitor program 
with regard to real-time operation. Specific references will be 
made to the SIGMA 5 computer, a product of Scientific Data Sys- 
tems, Inc., which has been proposed for this automation task. 


INTRODUCTION 


This report will attempt to demonstrate how the hardware and software 
capabilities of the SIGMA 5 computer can be utilized to develop a real-time 
control system for tracking station automation. When considering software 
capabilities, it should be emphasized that this is a preliminary report. The 
MONITOR program for the SIGMA 5 is currently undergoing product evaluation 
by SDS personnel who are independent of its developers. Consequently, SDS 
systems applications personnel do not yet know all the nuances of the BATCH 
MONITOR. However, from experience with other SDS monitor programs and an 
investigation of the preliminary BATCH MONITOR manual, the probable way 
that it will operate was determined. 


REAL-TIME SYSTEMS 

In order to explore the real-time capabilities of the SIGMA 5 computer, it 
is neccessary to define the characteristics of a real-time system. Basically, a 
real-time system contains an evaluation and control center (the computer) which 

*Miami University, Oxford, Ohio 
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responds to conditions existing within the system's environment. The control 
center performs processing based on these conditions such that commands to 
control the system are generated. To perform this task, several extensions of 
normal computer batch processing become necessary. 


Multiprogramming 

Since an environment is being controlled, it is necessary to have fast re- 
sponse to changes in the environment. This requirement places a great emphasis 
on the need for rapid input-output (I/O) to the computer. Further, there are 
situations when this requirement results in having more than one program share 
the computer. 

During a peak load, many messages may be seeking the use of the computer. 
The time between the arrival of messages may become so small that the process- 
ing time for the various messages may be greater than the inter-arrival time. 

It then becomes necessary to process the messages in parallel, i.e., have more 
than one program share the computer. 

Multiprogramming can be accomplished because the computer's functional 
units are able to operate independently and simultaneously. For example an 
input operation can be carried on simultaneously with a central processing unit 
(CPU) operation. If a program calls for data from an input file, the program may 
not be using the CPU while the information is being transferred to memory. 
Therefore, the CPU can be used by a second program. 


Interrupts 


While multiprogramming allows programs to share the computer through 
utilization of overlapping functions it becomes necessary in real-time systems 
for some messages to seize control of the computer regardless of which function 
the computer is performing. Thus, an interrupt control feature which seizes 
the computer and loads the correct program for message processing must be 
included in a real-time system. 

The following factors must be provided for effective interrupt control: 

A. Some priority measure must be attached to each message so that the 
computer processes the most urgent request. 
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B. The status of the program currently in control of the CPU must be 
saved when an interrupt occurs. 

C . The subprogram needed to process the interrupt must be identified and 
loaded from a file into high speed memory if it is not already present 
in memory. 

D. On occurrence of an interrupt of lower priority than the program cur- 
rently in control of the CPU, the interrupt message should join a queue 
for later processing. 

E. At the end of interrupt processing, control must be returned to the 
program that was current at the time of the interrupt. 

F. Interrupts should be able to be triggered (initiated) from several devices, 
e.g., analog to digital converters, input-output channels, the central 
processing unit, and timers. 


Memory Protection 

Having several programs which are being processed in paralled on the com- 
puter results in a problem of interaction which must be avoided. One program 
must be prevented from referring to memory locations used by another program. 
Such interaction would obviously impair the validity of either program's results. 

Further, if programs are to use a common subroutine, then the subroutine's 
logic should be such that intermediate results are saved when a program using 
the common routine is interrupted. When control is returned to the interrupted 
program, the intermediate results are restored and processing can continue in 
the common routine from the point in the logic at which the interrupt occurred. 
This type of programming makes the common routine re-entrant and is accom- 
plished by keeping the routine's instruction areas separate from the associated 
data area. When an interrupt occurs, the contents of the data area can be placed 
in temporary storage and the data area can be initialized for processing of the 
interrupting program. When the interrupt processing is complete, the temporary 
storage is loaded into the data area prior to resumption of processing on the 
interrupted program. 


Relocatability 


The scope of the system defines another characteristic of real-time 
computer applications. If all programs necessary for control can be resident 
in core at the same time, then, no provision for relocation of programs has to 
be made. In larger systems, programs may be called from file memory only 



when they are needed for processing certain messages. Since many combinations 
of programs may occupy core memory at different times, it is necessary for 
programs to be relocatable, i.e., different memory locations are used for 
storing the programs with each loading. Further, if memory space is at a 
minimum, it may become necessary to section a program such that the only 
logical portions of the program which are currently being executed are in mem- 
ory. The rest of the program is in library file storage and will be entered into 
core storage when needed for processing. 


Timers 


The most catastrophic occurrence in a real-time system is for a program 
to gain control of the computer and either go into a loop or experience a hard- 
ware failure such that the program never finishes its task. This results in 
locking out of all other processing tasks. To avoid this, the real-time system 
is provided with an internal timer which can be set so that no program can 
maintain control over the central processing unit in excess of the specified time. 

Another type of timer can be used to initiate interrupts at specified intervals. 
This interrupt could cause a sequence of programs to be initiated at these time 
intervals and thus, give the affect of cycling. 


IH. THE SIGMA 5 COMPUTER 

The SIGMA 5 computer is one manufacturer’s attempt at answering the 
requirements for a real-time system. Realizing the rapid response time 
necessary in real-time systems, computer hardware has been designed to meet 
these requirements. 


Input-Output 

Since a variety of messages may desire access to the computer concur- 
rently, a multiplexor input-output processor (s) (MIOP) is available which will 
handle up to 32 input-output devices simultaneously. A maximum 1-0 rate of 
500,000 bytes per second can be handled by each MIOP. The NTTF contract 
calls for two MIOPs. 
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A selector input-output processor (SIOP) is available but not specified in 
the NTTF contract. This processor is designed to quickly handle large amounts 
of data in a single burst. One input-output processor may be selected from as 
many as thirty-two devices. The transfer rate is approximately 4,700,000 bytes 
per second. This type of input-output is excellent for file processing but seems 
to have little application at the tracking station. 


Working Memory 

The core memory for the SIGMA 5 computer is organized in modules. 

Each module is capable of handling 4,8,12 or 16,000 words of memory. A maxi- 
mum of four modules can be attached to the system. The NTTF contract calls 
for two modules, each with 16,000 words of memory. 

Each memory module has six-way access. This allows the module to be 
directly attached to both the central processing unit (CPU) and the input-output 
processors (IOP). A particular module can only be accessed by one device at 
a time; however, two modules could be active concurrently by having one attached 
to the CPU and the other attached to an IOP. Thus, the computer could be per- 
forming an input or output operation using one memory module while simulta- 
neously performing calculations using information from the other module. 

By interleaving the addresses of memory words on the two modules, further 
simultaneity of processing can be achieved. Interleaving refers to the technique 
of assigned odd addresses in one memory module and even addresses in the 
other module. 

When a program is being executed, segments of the program logic may be 
executed in a sequential fashion, i.e., one instruction will be retrieved from an 
even location and executed and the next instruction will be retrieved from an 
odd location and executed, the next instruction from an even location, etc. 

The first consequence of interleaving is that, when the program is executing 
sequentially, the IOP would never be locked out of a memory module for a time 
greater than the cycle time.* This is true since all the even core addresses 
are on one memory module and all the odd core addresses are on the alternate 
memory module. Thus the CPU is accessing alternate memory modules when 
processing the program and the idle memory module is available for IOP 
selection. 


*Cycle time is the time to retrieve an instruction from core memory, bring it to the instruction 
register in the CPU, and rewrite the instruction in memory. Cycle time on the SIGMA 5 computer 
is .85 microseconds. 
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The execution of some programming instructions requires an alteration 
of core memory. Storing data from a register into core memory is an example 
of such an instruction. If the address of the memory location of where the data 
is to be stored is found in one memory module and the address of the next 
program instruction is in another module, interleaving increases processing 
speed since data can be stored while the next instruction is simultaneously 
retrieved and brought to the CPU. 


Real-Time Capabilities-Interrupts 

The SIGMA 5 computer uses interrupts to initiate its real-time capabilities. 
These are arranged such that when an interrupt is triggered, the computer 
automatically executes the instruction stored at a unique core memory location 
associated with the interrupt level. Thus, by storing a program control instruc- 
tion (such as an "Exchange Program Status Double-word" instruction - XPSD) 
at each of these interrupt locations, the routine to handle each particular inter- 
rupt can seize control of the CPU. 

It was observed when describing real-time systems that priorities must be 
assigned to messages. This is accomplished through hardware on the SIGMA 5. 
Each interrupt level has a unique core location associated with it and follows 
the rule; the lower the core location assigned to the interrupt level, the higher 
the priority that is associated with the interrupt. 

Interrupts can be classified as internal and external. Internal interrupts 
are triggered by such things as expiration of an internal interval timer, by a 
program request for input or output, etc. External interrupts are triggered by 
messages from devices on the MIOPs. Thus, interrupt levels over which the 
programmer has control can be grouped into clock, I-O, and external. On the 
SIGMA 5 computer specified in the NTTF contract, there are two levels of 
interrupts associated with the clock group, two levels associated with the 1-0 
group, and thirty-two levels associated with the external group. 

The interrupt level can exist in any of several conditions during the execu- 
tion of a real-time job. This allows assignment of priorities to vary somewhat 
even though they are tied to a fixed core location. Interrupt states can be defined 
using the following terms; 

1. A group of interrupts may be inhibited or active. 

2. An interrupt level may be enabled or disabled. 
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3. In addition to being either enabled or disabled, an interrupt may be in 
one of four mutually exclusive states; disarmed, armed, waiting, or 
active. 

When an interrupt level is enabled , the interrupt is waiting its turn for time 
on the CPU to execute the program associated with the interrupt. When it will 
be processed will depend upon its priority level. If an interrupt is disabled , it 
is removed from the priority sequence and must be enabled by a program in- 
struction before it can gain control of the computer. When a particular interrupt 
level is disabled a lower priority interrupt is permitted to gain control of the 
computer. 

When an interrupt level is disarmed and a signal is received on that level, 
no access to the computer is permitted and no record is kept of the signal that 
attempted to trigger the interrupt. This mode of operation could be used during 
testing of single programs so that there would be no interference with the test 
from outside sources. 

An interrupt level can be armed through program commands such that signals 
to that interrupt level will be accepted and the level will move to the waiting state. 


Once in the waiting state, the interrupt level must satisfy four conditions 
before the program associated with the interrupt level gains control of the CPU. 
These conditions are: 

1. The level must be enabled. 

2. The CPU must be at an interruptible point in its cycle. 

3. The group inhibit must be off. 

4. No higher priority interrupt level is in control of the CPU or is in the 
waiting state and enabled. 

Once the interrupt gains control of the CPU, the interrupt level is said to 
be active. The interrupt level remains in the active state until it is interrupted 
by a higher priority interrupt or until the program associated with it is completed 
and control is turned back to a supervisory routine. This can be done by having 
the program execute a ’’Load Program Status Word" (LPSD) or a "Write Direct" 
(WD) instruction at the logical end of the program. At this point the interrupt 
level reverts back to either armed or disarmed state. 

The disarmed, armed, waiting, and active states of an interrupt level are 
mutually exclusive. If an interrupt level is in waiting or active state and a 
signal to the interrupt level is effected it is lost. There is no queuing of 
requests for a given interrupt level. However, through the use of the waiting 
state, there is a queue of interrupts which are on different priority levels. 
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Real-Time Capabilities - Multi- Programming 


According to Scientific Data System’s personnel, there are no multi- 
programming capabilities as described in Section II of this report available on 
the SIGMA 5 computer. It is not possible for more than one program to share 
the same functional units of the computer. 


Real-Time Capabilities - Computer Control 

Because of external interrupts there is an overhead time for changing 
programs. SIGMA 5 hardware makes provision for rapid change of computer 
control and data areas. 

The control of the computer is centralized in a series of registers which 
comprise 64 bits of information. These registers are referenced collectively 
as the program status double word. The entire contents of this double word 
can be saved and reinitialized on five micro-seconds by a single program 
instruction. Theoretically then, control of the computer can be switched from 
one program to another if five micro-seconds. 

The registers in the status double word that are of interest for the purposes 
of this report are: 


1) IA Introduction address 
RP Register pointer 

CI 1 

II > Interrupt group inhibits 

EI J 

MS Master/slave mode 
WK Write key 


The IA is a seventeen bit register which contains the address of the next 
program instruction to be executed. When a program is being executed, this 
register is increased by one each time a program instruction is retrieved from 
storage. Thus, the instructions are executed in sequence except when a branch 
in the program logic is encountered. When this takes place, the core location 
to which the program is to branch is transferred to the IA. 

Since the I A register controls which instructions are executed, this is the 
register which must be changed when a new program is to be executed due to 
an interrupt. This change is accomplished by having an ’’Exchange Program 
Status Doubleword" (XPSD) instruction located at the core location associated 
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with each interrupt level. The operand of the XPSD instruction contains the 
location of a new program status doubleword which is associated with the 
interrupt. This program status doubleword is exchanged with the current program 
status doubleword. Thus, the I A changed and the first instruction of the program 
associated with the interrupt is executed. 

In addition to control being switched, some provision must be made for 
protecting the data associated with the old program and for initializing the data 
areas for the new program. To a degree, this is also accomplished through the 
program status doubleword. 

There are groups of sixteen general purpose registers contained in the 
CPU. The NTTF contract calls for two such groups but up to sixteen such groups 
are available on the SIGMA 5. These general purpose registers can be used as 
accumulators, index registers (seven of the sixteen registers), temporary storage, 
control information, counts, pointers, etc. This information is relevant to a 
particular program and must be saved when the program is interrupted. 

Within the program status doubleword, the register pointer (RP) is used to 
point to the block of sixteen registers currently in use. When the program status 
doubleword is loaded by the XPSD instruction, a new value for the RP can be 
loaded which will point to a new register block. Thus, volatile data associated 
with the interrupted program doesn't have to be stored but can be preserved by 
merely having the register pointer changed. Since the execution time for the 
XPSD instruction is five microseconds, the entire enrivonment of the computer 
can theoretically be changed in this amount of time. 

In reality, with the NTTF system, this five microsecond switching will not 
be possible if the MONITOR system is used. The MONITOR requires a register 
block for its own use and therefore, only one register block is available for 
application programs. Thus, if the entire register block is used for each appli- 
cation program, the register block will have to be stored each time one program 
is interrupted and another program is loaded. The register block will then have 
to be initialized for the new program. 

The storage and initialization of the entire register block could be accom- 
plished by the first routine in the program associated with the interrupt. A 
"push down" instruction can store the contents of the registers in a group of 
storage locations reserved for this purpose. Thus a queue of the contents of 
the registers for interrupted programs is formed. A "load register" instruction 
can be executed so that the contents of an initialization area are loaded into the 
general purpose registers. The complete switching of all sixteen registers would 
take approximately 36.5 microseconds. 
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Protection of data areas that are not in the general purpose registers but 
are in core memory can be accomplished by the use of memory wire protect. 
This allows memory to be protected in 512 word pages. A two bit hardware 
"lock" in the CPU is associated with each of these pages. The program status 
double word contains a two bit write key (WK) which is checked against the 
lock before core memory can be altered by a program instruction. If a violation 
of the write protection occurs, the memory is not altered and control of the 
computer is sent to location X’40’ where a branch to a programmed procedure 
to handle this error should start. 


IV. MONITOR PROGRAMS 

As an introduction to the real-time capabilities of the SIGMA 5 computer, 
the characteristics of real-time systems were described. It will help with the 
exposition of the BATCH MONITOR if monitor programs are also explained in 
general. 

The monitor program for a computer is provided as a controlling program 
for smooth hardware operation. The program normally consists of five major 
sections which are the supervisor, translators, the load routine, the foreman, 
and the input-output control system. 


Supervisor 

At least a portion of the supervisor program is in core memory at all 
times. This section of the monitor program provides for communication 
between the programmer and the computer system and for recovery from errors 
during program execution. 

Communication between the programmer and the computer is accomplished 
by a series of control cards. Each program can be identified and the monitor 
program processing to be performed on the program can be defined by these 
control cards. Key words which the supervisor program recognizes initiate 
a chain of events. Examples of these events are: 

X. Dynamic assignment of input-output devices to the files used by the 
application program; 

2. Calling of the assembler or compiler programs into memory; 

3. Calling the loader routine into memory; 

4. Keeping accounting data on program translation and execution times. 
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Uninterrupted processing can be achieved since the supervisor handles 
error conditions. 

The application programs are stacked together as jobs to be done by the 
computer. When one program’s processing is completed, control is returned 
to the supervisor program so that the next program job can be called and started. 
Should an abnormal or error ending occur within the application program, a 
branch to the monitor supervisor routine is attempted. If the supervisor program 
has not been destroyed as a result of the error, the supervisor program is able 
to recover and move on to the next program job in the sequence. 

This arrangement of the job stream of programs to be executed has given 
rise to the term ’’batch shop" processing, i.e., a batch of independent programs 
are grouped together for processing and the monitor's supervisor routine 
handles the housekeeping necessary for processing of these programs. It is an 
obvious extension of this processing for a scheduler routine to be programmed 
as part of supervisor routine. This enables the program jobs to be run to be 
assigned a priority level and the scheduler routine will cause the highest level 
priority programs to be executed first. 

The scheduler should not be confused with the priority interrupt levels as 
described for the SIGMA 5 computer. The interrupt priority levels are asso- 
ciated with hardware and real-time applications. The scheduler is associated 
with software and batch operation. 


Translators 


These programs are used to translate source language programs into object 
or machine code. Assembly programs are characterized by the fact that they 
generally translate one executable source program statement into one machine 
language instruction. Any computer of a significant size would have an assembly 
program associated with it so that programmers would not have to program the 
computer in machine language. The assembler program for the SIGMA 5 
computer is called SYMBOL. 

Compiler programs are characterized by the fact that they translate one 
executable source program statement into several machine language instructions. 
A compiler program attempts to standardize a source program language for 
several computers. A program written in a compiler language could be trans- 
lated (with minor changes) into the machine language for several different 
computers. Examples of compiler type translators are FORTRAN, COBOL, 
ALGOL, and PL-1. 
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The translator programs are part of the monitor system. They would be 
stored on some auxiliary mass storage device such as disk storage and called 
into high-speed core memory only when they are needed to convert a source 
program into machine language. 


Loader 

r 

The loader routine is called into high-speed core memory from a mass 
storage device by the supervisor program when it is necessary to load machine 
language programs into high-speed memory. This could be done immediately 
after translation of the source program into machine language. This is the 
"translate and go" mode of operation and is excellent for testing program logic. 
If errors in logic are found in the program, corrections can be made to the 
program in source language. The program is then translated, loaded, and run 
with test data again. The "translate and go" mode of operation does require a 
rapidly executing translator so that all the computer time is not wasted on 
translations. 

The loader routine will also load machine language programs which have 
been generated previously and punched on cards or stored, in auxiliary mass 
storage. This would be the normal mode of operation when the programs are 
operational, i.e., the debugging phase has been completed. Having the programs 
in object form eliminates the time-consuming translation from source language 
to machine language each time the program is executed. 

A major function of the loader is to bring all subroutines and library 
routines which were requested by the main application program into memory 
and link them to the application program. The loader will also perform the 
relocation of the routines brought into memory. 


Foreman 


This monitor routine is responsible for the correct functioning of the 
computer during application programs for either real-time or batch processing 
applications. Obviously then, it must be resident in core memory at the time 
of execution of the programs. 

In batch processing and real-time applications, this routine is responsible 
for the operation of the input-output devices. It gives the actual command that 
initiate movements of the input-output mechanisms. Further, the routine 
monitors the performance of device. If error occurs during the input-output 
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operation (e.g., a parity error when reading data), the FOREMAN routine would 
cause the 1-0 device to perform the operation again. 

Occasionally programs will be "trapped" during execution. This can result 
from several situations. An example would be a trap resulting from a division 
by zero. A trap could also be caused by a floating point data word which over- 
flows its exponent range. When traps occur, it is the responsibility of the 
FOREMAN routine to perform the appropriate logic associated with the trap 
condition. 

In real-time applications, an addition function is demanded of the FOREMAN 
routines. The routines must also handle external and internal interrupts. The 
interrupt service routines and the linkage to programs connected to the interrupt 
level are considered part of the FOREMAN. 


IOCS 

These routines are common to all input-output requests. To alleviate the 
tedious task of writing these routines each time a reference is made to a file 
of input or output data, these routines are supplied in the MONITOR package 
and can be referenced by a single source statement in the application program. 
These routines interact extensively with the FOREMAN routines for 1-0 control. 
They become an integral part of the program and thus, they are resident in core 
during the execution of the application program. 


V. SDS BATCH MONITOR 

The BATCH MONITOR for the SIGMA 5 computer will be available in 
October, 1967. This monitor will have a 4,000 word resident core load. The 
supervisor, translator, IOCS, and loader sections will be provided. The programs 
necessary for the real-time FOREMAN control will have to be written by the 
user. These real-time capabilities are termed to be FOREGROUND programs. 
Batch processing can be accomplished when the SIGMA 5 is not in the real-time 
mode. This batch processing is termed to be BACKGROUND programming. 

It will have to be determined whether the foreground programs to handle 
various equipment interrupts will need to be resident in core or whether they 
can be stored in disk storage and called into core only when needed due to a 
particular interrupt. This decision will depend upon the quantity and size of 
programs necessary to perform all requirements for the NTTF station, the 
size of computer memory, and the speed of response to various interrupts. 
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Further, a basic decision must be made concerning the language in which 
the foreground tasks are to be programmed. This will have an effect on storage 
and timing requirements. The remainder of the report will discuss the advan- 
tages and disadvantages of the METASYMBOL and FORTRAN languages and how 
the interrupt routines are made part of the FOREMAN routine in the monitor 
package. 


Meta-Symbol 

META-SYMBOL is the assembler program that will be available with the 
BATCH MONITOR in October, 1967. This language gains much of its power 
from the fact that user defined macro instructions are easily generated during 
an assembly. By the use of a "procedure directive" instruction (PROC), user 
defined macro instructions can be declared for use later in the executable 
program. The user can define a general routine to perform some processing 
function which is recurrent within his program. When the particular processing 
function is needed in the executable program, the name of the procedure directive 
is given in the command field (in the same way as a source language operation 
code would be given) and the field elements to be used in the procedure are 
supplied as arguments. This source coding causes machine language instructions 
for the routine to be generated in line with the machine language for the rest of 
the program. The generated machine code will have the correct addresses for 
the specified field elements. 

META-SYMBOL also has the advantage of being able to reference monitor 
subroutines which are supplied with the BATCH MONITOR. These routines are 
routines which are helpful to the real-time systems programmer. They include 
routines that: 1) can set interval timers so that a fixed cycle of operations may 
be set in progress due to a timer interrupt: 2) provide for dynamic allocation 
of interrupt levels: 3) control reentrance to subroutines by moving context data 
from the subroutine to a dynamic storage area and returning the data to the 
subroutine after the higher priority interrupt has been serviced. These features 
and others are supplied with the BATCH MONITOR and are resident in core. 

A reference to them by a CAL instruction within the source program causes a 
branch to the monitor routines to be generated. 

By choosing to program the foreground routines in META-SYMBOL, the 
user would gain the advantage of being able to begin writing and debugging 
routines immediately. The advantages of procedure directives would enable 
good systems programmers to provide a powerful assembler. MONITOR 
subroutines offer a convenience in dealing with real-time system demands but 
should be investigated in detail to see if they are efficient enough for the time 
requirements of the NTTF tracking system. 
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The disadvantage of META-SYMBOL as a programming language for the 
tracking programs is the fact that it ties the programs to the SIGMA 5 computer. 
Should another computer be specified for other tracking sites, then the programs 
would have to be re-written. For this reason, the advantages of compiler level 
programming for the foreground programs must be investigated. 


FORTRAN 

The delivery shedule for the FORTRAN compiler leaves something to be 
desired when a real-time application is being considered. FORTRAN is 
scheduled to be delivered in four versions beginning in November, 1967 a nd 
ending in January, 1969. The special features needed to make FORTRAN easily 
connect with a real-time operating system will not be available until the January, 
1969, version. The most glaring omission from the earlier versions will be 
the lack of a reentrant feature for subroutines. The monitor's capabilities 
cannot be utilized for this feature since no communication between the generated 
user program and the BATCH MONITOR routines is supplied as it is with the 
META-SYMBOL assembler. 

Of course, the user can write his own routines so that any subroutine can 
become re-entrant. Upon being interrupted, it would be the responsibility of 
the FOREMAN program to save the program status double word, the current 
general purpose registers, the contents of location X(4F) and any volatile data 
that might be affected by the program servicing the interrupt. The last require- 
ment would refer to any common user written subroutines that service interrupts, 
and FORTRAN library programs such as SINE, COSINE, etc. Since earlier 
versions of FORTRAN will not have separate areas for data and instructions in 
the library subroutines, it may be necessary to duplicate library routines and 
include some type of software pointer in the user written FOREMAN routine so 
that an available library routine for a particular function can be assigned to each 
unique interrupt message. 

The obvious advantage of a compiler language over an assembler language 
is the ease of writing and debugging the program. In this case, some of this 
advantage is lost because of the delivery dates of the real-time FORTRAN 
compiler. On the other hand, more compact code is generated when good systems 
programmers use an assembler language as opposed to a compiler language. 

A study must be made to determine whether the real-time control program's 
memory requirements are loose enough to generate code by FORTRAN pro- 
gramming. Further, the scope and consequences of the project must be deter- 
mined. If there is a high probability that the final real-time control package 
will find its way into all STADAN stations, then it would be convenient to have 
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the programs in FORTRAN so that the control computer can be rebid. Finally, 
FORTRAN coding offers the advantage of being easily changed when the need 
arises due to various tracking experiments that will be developed in a testing 
and training facility. 


System Generation 


Regardles of which language is used to generate the foreground tasks 
associated with real-time control, they must be added to the delivered BATCH 
MONITOR system becoming a permanent part of the resident core load and/or 
an addition to the portion of the monitor that resides on the disk during opera- 
tion. These new routines would be added to the monitor during a systems 
generation computer run. 

It will be recalled from the section on monitors in general that the pro- 
grammer can communicate with the monitor through a series of control cards. 
Thus the monitor is able to expand itself by recognizing the control cards and 
adding the programs that follow to itself. 

The procedure to expand the monitor by attaching routines to interrupt 
locations would be as follows. The routines themselves would be in a source 
language (META-SYMBOL or FORTRAN) . These would be translated to machine 
language by the appropriate assembler or compiler program. The machine 
language programs would be stored on an element output (EO) file such that 
each had a unique name and each was a unique object module. 

The load routine which was described in the above section would now be 
used to create a new monitor by loading each object module into memory and 
attaching it to an interrupt location. 

Specifically, the ! LOAD control would be used. The element output file 
which was created during the translation phase above would now be referred to 
as the element input (El) file by the appropriate parameter on the ! LOAD card. 
Each object module that was to be loaded from the El file would be specified in 
the parameter argument list. Other files such as the System Library (SYSLIB) 
could be searched to have certain programs contained thereon added to the 
new monitor. 

Another control card(s) following the ! LOAD card would attach the object 
modules to a given interrupt location. The : INTR card specifies the name of 
the object module and the hexadecimal address of the interrupt with which it 
is to be associated. The function of this card is to load the XPSD instruction 
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in the correct memory location for the interrupt level. Thus, during program 
execution, an interrupt at this level will send control to the object program to 
handle the interrupt. 


VI. CONCLUSION 

It is difficult to determine exactly how each instruction or option within 
the BATCH MONITOR control language will operate or what its effect will be 
when used in conjunction with other control instructions without having the 
benefit of a computer system to test combinations. Hopefully, this report does 
clarify some of the underlying principles behind real-time control systems and 
that these principles will serve as a guide when actual writing and testing of 
the necessary programs becomes a reality. 
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ON THE USE OF ADAPTIVE TECHNIQUES TO IMPROVE 
ANTENNA SERVO ACCURACY 

by 

H. J. Perils* 


SUMMARY 

The basic antenna tracking system is analyzed under various 
conditions. The study suggests some desirable relationships to 
maintain within the framework of the existing configuration. In 
order to improve performance under different tracking environ- 
ments an executive form of adaption is recommended. This em- 
ploys available external information in a scheme to generate an 
adjustment vector in the tracking servo system. 


INTRODUCTION 

Under the sponsorship of the 1967 ASEE-NASA Faculty Fellowship program, 
the author initiated a study to improve the tracking accuracy of existing servo 
systems. Since the systems considered are currently operational and various 
large antenna systems have been operating for over there is a limited 

amount of design and operating data available. 1, 2,3 Hence, the work started 
with a review of the readily available design and operating manuals and reports. 
This was combined with visits and talks with tracking station personnel of the 
STADAN-Rosman I and the MSFN-Goddard NTTF facilities. 

Although it is claimed that sophisticated simulation and analysis schemes 
are available, quite reasonable analytical results are obtained with the use of 
a mixture of simplified linear lumped models, Bode diagrams, and some "rule- 
of -thumb” folklore. This is made possible by a combination of state-of-the-art 
specifications, careful system alignment and calibration, the large masses in 
the system, and the use of highly accurate sensors and powerful hydraulic servo 
motors. The delivery schedule, the low-pass nature of the system, and the lack 
of reliable numerical values for stochastic parameters makes the employment 
of elementary design and error analysis methods understandable. 


*Rutgers, The State University, New Brunswick, N. J. 
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The servo system design in auto-track mode was reviewed. The nature of 
the sources of auto-track error was analyzed. The objective of this work was 
the development of ideas for an executive-type adaptive adjustment schedule. 
Within the basic framework of a fixed configuration a set of adjustable parameters 
was selected for each axis and a set of information parameters was selected. 

The information parameter vector is employed in an adjustment computer to 

generate the adjustment parameter vector. 

/ 


THE SYSTEM ANALYSIS OF THE TRACKING SERVO 

The antenna tracking servo system studied was of the X-Y type, specifically 
the MSFN 30-foot system. 4 To simplify the analysis, the following discussion 
is addressed to a single channel. Since the channels operate independently and 
are of similar composition, the X channel is discussed because it has the greatest 
load. 

At low frequencies, one can treat the load as a pure inertia J (reflected back 
to the motor shaft), and the hydraulic servo actuator subsystem can be represented 
by a torque constant K m . A linear, low-frequency model of the tracking servo is 
given by Figure 1. The transfer functions G c (s) and G Cr (s) are the position 
and rate operational amplifiers with their associated networks. The motor-to- 
antenna gear reduction is given by N g , and the tachometer gain is shown as K t . 

The input closed-loop position command is ( 9 i . It includes both a desired angle 
signal 6 X and a contaminating noise signal n. The system is subject to a wind 
disturbance. This is introduced as a torque signal T w which reflects back to the 
motor shaft. 

The output antenna position 6 k is expressed in terms of the responses of 
all of the inputs. Thus, in Laplace transform notation 


let 


S A ( s > 






N (s) + 



(s) 


( 1 ) 


F t (s)-^-(s) 
i 

9 a 

F 2 (s)--i-(s) 
1 w 


( 2 ) 

(3) 


E(s) = S x (s)-6? a (s) 


(4) 
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Then, 


E (s) = (1 - Fj (s» .0, (s) - F ;1 (s) N (s) - F 2 (s) T w (s) (5) 

At this point, one might consider some linear optimization schemes. If a 
semi -free configuration is assumed along with stochastic signals for 9 X , n and 
T w , a Wiener filter can be< found to minimize the RMS value of the error. 5 A 
more realistic approach would be to recognize the fact that 6 X is deterministic 
and can be approximated by a time function. 6 Thus, the extension to the Wiener 
filter for a restricted class of mixed inputs could be used. 7 

Since the system being analyzed actually exists, it would seem more approp- 
riate to consider a fixed configuration studied under some simplifying conditions. 
Thus, the following two situations will be separately analyzed: 

(a) 9 X large 

(b) 9 x small 

Case (a) could represent a low altitude orbiting space vehicle with a large 
9 . In such a case the error in Equation (5) would be dominated by the first 
term. That is, the system error e (t) could be approximated by the actuating 
error e a (t) , and a classical deterministic analysis would follow, 8 

Case (b) could represent a high altitude, slowly moving space vehicle. In 
this case the noise and wind disturbance terms would dominate Equation (5). The 
problem could be treated as a fixed configuration Phillips’ -type optimization in- 
volving some adjustable parameters. 9 

In any event one must work with Figure 1 and develop the expressions for 
F x (s ) and F 2 (s). 

Thus, 


where 


(s) = 


G(s) 

1 + G(s) 


( 6 ) 


G (s) 


J-G c (s) G (s) 
Ng c p 

s 


( 7 ) 
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and 


G r (s) 


K 

fv s) 

s +K t 5iG Cr (s) 


( 8 ) 


Since a low frequency study is being employed which ignores the high fre- 
quency dynamics of the antenna and the hydraulic actuator, the G c (s) also can 
be treated as a constant. That is, r 


Let 


K 


r 



Thus, 


G r (s) 



s +K t K 


r 


In the auto-track mode the G c (s) is of the form 

c p 


( 9 ) 


( 10 ) 


G c (s) 

p 


K p ( s +"1) 
S 


( 11 ) 


Thus, the over-all system reduces to a simple unity negative feedback Type 
2 system. The open-loop transfer function is given by 


G (s) = 


^ k , k p< s+ “i> 
s 2 (s+K t K r ) 


( 12 ) 


or in time -constant form 


K (1 + s/^> 

G (s) =— — 

s 2 (1 + s/w 2 ) 


(13) 
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where 


<V=K t K 


(14) 


K ■ lim 

; a 

s— * 0 


[s 2 G (s)] = 


Ng K t 


Cl>, 


(15) 


For s = j co one can sketch the Bode asymptotic amplitude response of G ( j co ) . 


v K (1 + j co/o) t ) 
G(jw) = 

(j") 2 (1 + j "A> 2 ) 


(16) 


Figure 2 is a representation of the Bode asymptotic amplitude response for 
a typical stable system of this form. In general, the lead network corner is 
greater than 1 radian/sec., and for stability the co 2 given by Equation (14) must 
be larger than co. . For a reasonable relative stability co c is normally located 
between co t and co 2 . 

Simple geometry will show that 

K = co co. (17) 

a c 1 

Also, for maximum Phase Margin it is easy to show that the crossover fre- 
quency co c should be at the geometric mean of the two corner frequencies. That 
is, 


CO = V CO. CO. 

C 1 z 

Let 

"a 

k= — , > 1 

Then, by combining Equations (13), (17), and (18) 

o? (s + ct) /i/k) 

G ( s ) = — — 

s 2 (s + vie a> c ) 


(18) 


(19) 


( 20 ) 
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If a normalized s -plane is defined as the A -plane and 


a - /k 

( 21 ) 

s/co = k 

C 

( 22 ) 

a 

G (\) - V a ' 
k 2 (k + a) 

(23) 


For such a system by standard geometric and trigonometric relationships 
one can show the constellation of poles and zeros in the X -plane by Figure 3. 
From this diagram it follows that 


£ = 


a - 1 
2 


(24) 


v = - £ 2 


(25) 


Since the closed-loop has a real pole at A. - - 1 , the system cannot be treated 
as a simple second order system or even a dominant pole system. That is, the 
dipole separation must be considered. A fractional closed-loop dipole separation 
taken with respect to the pole at A = - 1 is given by 

DS=lzi (26) 

a 

In the real frequency domain the Phase Margin can be written as 

PM = 2 tan" 1 a - 90° ( 27 ) 

As k increases the Phase Margin and the damping ratio £ increase, and 
the dipole separation and the normalized damped natural frequency 17 decrease. 
These variations are given in Table 1 for the interval 1 1 a 5 3. 
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Figure 3-The Open and Closed-loop Poles and Zeros in the X— Plane 


Table 1 


k 

a 

PM(°) 

i 

V 

DS(%) 

1 

i. 

0. 

0. 

1. 

0. 

2 

1.414 

19.4 

0.207 

0.976 

29.2 

3 

1.732 

30. 

0.366 

0.930 

42.1 

4 

2. 

36.8 

0.5 

0.865 

50. 

5 

2.236 

41.8 

0.618 

0.785 

56.4 

6 

2.449 

45.4 

0.725 

0.686 

58.0 

7 

2.646 

48.6 

0.823 

0.565 

62.2 

8 

2.828 

51.0 

0.914 

0.405 

64.5 

9 

3. 

53.0 

1. 

0. 

66.7 


The unit step response of this third-order system can be found via standard 
inversion tables. 10 Using normalized time r, the response over the open interval 
1 < a < 3 is given by 


e ( T ) = 1 + ( a ~ 1 ) 
aW 2 ( 1-0 


— T ^ 

e T + — 


±-!£ + 1 

a 2 a 


1/2 


2 ( 1-0 


'^ T sin ( ’j)T + a) 


T > 0 

1 < a < 3 


where 


a = tan 1 —H — - tan -1 — - tan -1 -2- 

I- £ 1-4 -4 

a 


(28) 


r = co t (29) 

C 

The response given by Equation (28) was computed for k = 1.5, 2.0, 2.5, . . . , 
8.5 using a digital computer at Rutgers University. Figure 4 is a plot of the re- 
sponse family for k = 2, 4, 6, 8. The effect of the dipole separation is quite 
apparent. For low damping ratio the dipole separation is small, and the system 
resembles the simple second order case. Increasing the damping ratio lowers 
the peak overshoot. The effect of the zero is a tendency to increase rise time. 
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Figure 4— A Family of Unit S 


The effect of the dipole is to "smear out" the backswing. Since the gain, dipole 
position, etc. are all dependent upon k, the system has the same "response time" 
(time to the first peak) t which is approximately 

t (30) 

P 00 
c 

< 

This is less than the equivalent second-order case because of the zero's effect 
on the leading edge response. The dipole effect on the backswing can be seen in 
Figure 4. It is interesting to note the approximate common crossing points at 
co c t ~ 1.5 and 4.5 which are the results of the presence of the dipole. Table 2 
gives a more complete listing of peak overshoot PO and normalized (±5%) settling 
time a> c t s (5%). The backswing smearing causes the settling time to increase 
with k for k > 4.5. 


Table 2 


k 

PO(%) 

% t s (±5%) 

1.5 

80.0 

25.8 

2.0 

67.6 

13.8 

2.5 

58.8 

10.7 

3.0 

52.3 

7.7 

3.5 

47.2 

7.8 

4.0 

43.2 

7.4 

4.5 

39.9 

5.3 

5.0 

37.1 

5.4 

5.5 

34.8 

5.6 

6.0 

32.8 

5.7 

6.5 

31.0 

5.8 

7.0 

29.5 

6.0 

7.5 

28.2 

6.1 

8.0 

27.0 

6.3 

8.5 

25.9 

6.5, 


Throughout the above work, it is seen that, under the maximum Phase 
Margin constraint, the system is completely described by k and a> c . From the 
point-of-view of maximum steady-state position error, one can write 


K 


( 31 ) 
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or 


K 


6 2 

x CO 

max. ^ c 


min. 0 


V^k 


(32) 


Thus, to minimize static error one requires a high value of low frequency 
gain. This is reflected in expression of Equation (31). This will be the dominant 
consideration in the Case (a) situation. 

In the Case (b) situation E (s ) ~ E x ( s) where 

E x (s) = - (F 1 (s) N (s) + F 2 (s) T w (s» < 33 > 

F x (s) is the same as before and is given by Equation (6). F 2 (s ) is related 
to Fj ( s ) by the following 


F 2 ( s) = Ms) 


(34) 


where 


G (s) = N K G (s) G (s) 

c v ' g mc v/ c v/ 


(35) 


Both n (t ) andT w (t) are ergodic stochastic signals, and, if the following assump- 
tions are made, 


<n (t)> = <T W (t)> = <n (t) T w (t +r)> = 0 


(36) 


the spectral density expression is given by 


< ‘.,. 1 ( s )= F . <-s)F 2 (s) 


$ (s) + 

n n x / 


^ww ( s ) 


G r (-s) G (s) 


(37) 


Using the low frequency approximation of Equation (9), 


G c ( S ) = K W M 1 


(1 + s/tUj) 


(38) 
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whore 


K w = N J K K 

W g J r p 

1'orG (s ) written in the following form 

K (1+s/^) 

G (s ) = - 

s 2 (1 + s/o> 2 ) 

it follows that 


(39) 


( 40 ) 


K.d+s/^) 

( S ) 

s 3 /co n + s 2 + K /oo. s + K 

2 . ala 


F a (s') _ K w Wj 

G C < s ) s 3 / co„ + s 2 + K ; u s + K 

£ 3 1 3 


( 41 ) 


( 42 ) 


where 


K 


k W N 2 K t K r J 


1 K 

K =__E.w 1 
a N K. 1 

g t 


( 43 ) 


( 44 ) 


Therefore , 


E, (s) = 


K a; n 


- 1 


K 


+ 1 


(1 + s/coj) N(s) + 


K w 


T w (s) 


( 45 ) 


CO, 


The noise can be considered to be approximately "white”. The literature of 
wind torque disturbance on large antenna makes use of the following approxima- 
tion (at least as a first order) for the spectral density 11 ’ 12 
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(46) 


«w oo 


(1 - s 2 /p 2 ) 


(lb -ft) 2 
rad. /sec. 


Where is proportional to the fourth power of the wind velocity. Using this 
and the following form for the noise 


B" 


4> (s) = 

e i e i v ' 


0 (s) 

nn v 7 

= B 2 >0, 

(rad. ) 2 

rad. /sec. 

fs- / 

— + — 1 

' B S ' 

CN 

i 

v 2 cJ 2 

< V 

- B 2 K S. 


(47) 


+ 1 


s 4 + _1_ l + ± 

+K a +W 2 


3/1 1 

s + k + - 


CO^ V 


s 2 + (— +.— ] s + 1 

\V 0), 


(48) 


[*] 


Where [ * ] denotes the property that this factor is identical to the other 
denominator factor with s replaced by (- s). 

By Parseval’s Theorem 


< e i (!)/> = 


/•J 00 

_L f $ 

J_ jm ^ 


(s) ds 


(49) 


Substituting Equation (48) into (49) results in a standard Phillips’ type inte- 
gral with n = 4. These integrals are tabulated. 5 The resultant expression can 
be written in terms of a set of normalized coefficients. 


2 <e 2 (t ))> 
v B 2 

n 


L + /3 + l 

i-/jk §___ 

1 +k/3 ‘ 


1 +P 


1 

+— 


1+/3 l+k/3 

l+k/3 %k(l+/3) 




B 2 > 0 

n 


(50) 



Where 


-\ 

k = o> 2 / co % 
fj. = ujv 


P - o> x / v 


> 


y 



(51) 


To obtain a feel for the variations of Equation (50), digital computer runs 
were made at Rutgers for various sets of normalized coefficients. The maxi- 
mum Phase Margin approximation was not made in the above work. However, 
to make sense, p should be in the range fi/k < p < /x. A convenient value for 
P is the geometric mean of the two limits. That is, 



(52) 


This is identical to the maximum Phase Margin condition. 

The results of digital computer rims are tabulated in Table 3 for sets of y 
and /x values against k , using the maximum Phase Margin P . The following 
trends are noted: 

(a) for all cases a minimum exists at k = 4 

(b) with no wind the RMS error is dependent of /x 

(c) with wind the RMS error decreases with a decrease in y 

(d) with wind the RMS error goes down with an increase in /x 

(e) for a given /x and k the error decreases with a decrease in y. 

These normalized quantities called F in the table must be multiplied by a 
factor. That is, 


<e\ (t)> “ y R2 F (k, fx, P, y) 


(53) 
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Table 3 


r(y s.L) 


r(y = 1 . 5 ) 


X 

0.25 

2 

4.83 

4 

4.00 

6 

4.14 

8 

4.38 

10 

4.62 


T(y = 2. ) 


X 

0.25 

0.5 

1 

2 

4 

2 

7.46 

7.04 

6.24 

5.45 

5.03 

4 

6.44 

6.03 

5.33 

4.63 

4.23 

6 

6.80 

6.35 

5.61 

4.87 

4.42 

8 

7.28 

6.78 

5.99 

5.20 

4.70 

10 

7.76 

7.22 

6.38 

5.54 

5.00 


X 

0.25 

0.5 

1 

2 

4 

2 

6.14 

5.93 

5.54 

5.14 

4.93 

4 

5.22 

5.02 

4.67 

4.32 

4.11 

6 

5.47 

5.24 

4.87 

4.51 

4.28 

8 

5.82 

5.58 

5.18 

4.79 

4.54 

10 

6.19 

5.92 

5.50 

5.08 

4.81 


T(y = 

: 2 . 5 ) 


X 

0.25 

0.5 

1 

2 

4 

2 

8.77 

8.14 

6.95 

5.76 

5.13 

4 

7.66 

7.05 

6.00 

4.95 

4.34 

6 

8.13 

7.45 

6.34 

5.24 

4.56 

8 

8.73 

7.98 

6.80 

5.62 

4.87 

10 

9.33 

8.52 

7.26 

6.00 

5.19 


T(y = 3 . ) 


T(y = 3 . 5 ) 



0.25 

0.5 

1 

2 

4 

2 

10.08 

9.25 

7.66 

6.07 

5.23 

4 

8.88 

8.06 

6.67 

5.27 

4.46 

6 

9.46 

8.55 

7.08 

5.60 

4.70 

8 

10.19 

9.18 

7.61 

6.03 

5.03 

10 

10.90 

9.81 

8.14 

6.46 

5.37 


X 

0.25 

0.5 

■ 

2 

4 

2 

11.40 

10.35 

8.36 

6.38 

5.33 

4 

10.10 

9.08 

7.33 

5.59 

4.57 

6 

10.79 

9.66 

7.81 

5.97 

4.84 

8 

11.64 

10.39 

8.42 

6.45 ' 

5.19 

10 

12.47 

11.11 

9.02 

6.92 

5.56 


T(y = 4 . ) 


X 

0.25 

0.5 

2 

12.71 

11.46 

4 

11.31 

10.10 

6 

12.12 

10.76 

8 

13.09 

11.59 

10 

14.04 

12.41 



2 

4 

6.69 

5.43 

5.90 

4.69 

6.33 

4.97 

6.86 

5.36 

7.38 

5.75 


L39 






































It has been shown experimentally that v is a constant and is independent of 
wind velocity. 11 Typical values used are v = 0.11 radians/second 11 and v = 1.0 
radians/ second. 1 2 


If v= 1.0 rad./ sec. is assumed, a possible desirable set of nominal condi- 
tions might be : 


-\ 


k % 4 
co > 2.0 rad. /sec. 

C 

PM ~ 35° 


> 


(54) 


y < 2.0 




For this set 


/ B a \ 1/2 

<56> 

The effective wind stiffness constant (lb-ft/rad) K w is given by Equation (39). 
If one assumes that , K p , K t , and K r are adjustable, it should be possible to 
obtain a compromise for each set of B* and in the Case (b) situation. Simi- 
larly, if the situation is of the Case (a) type an appropriate K a given by Equation 
(15) can be found to satisfy Equations (17) and (32). 


THE USE OF EXECUTIVE CONTROL 

Based upon the work of the previous sections, the notion of self adaption is 
brought to mind. Since the tracking servo does not have to be self-contained, 
available outside information can be used to help to adjust the system in accord- 
ance with the objectives set up in the previous section. Such a system could use 
the concepts of executive control. 13 That is, both measured values and predic- 
tion data can be used to develop the needed Case (a)/ Case (b) decision and the 
associated set of required adjustments. 

Typical predictions available by computation or from past experience are: 

1. Minimum expected range 

2. Maximum expected Y 
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3. Maximum expected X 

4. Expected AGC level 

Typical measurements available for up-dating: 

1 . RMS value of wind velocity 

2. Average value of AGC 

The predictions, measurements, together with nominal relationships and 
values for the system parameters can be used to develop the required adjustment 
vector. A conceptual block diagram of the executive-control adaptive tracking 
servo system is given in Figure 5. In the limited time available to devote to this 
project it was not possible to specify the complete details for the adjustment 
computer. However, it is felt that the proper directions and much of the basic 
system analysis has been performed. Thus, this work can be used as a start 
toward the realization of antenna tracking servo self-adaption as part of an over- 
all station automation program. 



Figure 5— The Over-all Executive Control Conceptual Block Diagram 
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CONCLUSION 


The overall system of an antenna tracking servo was analyzed in detail. The 
study was divided into two parts: (a) the case in which requirements for antenna 
acceleration dominate, and (b) the case in which wind gust loading and tracking 
receiver noise dominate. The study indicated the desirability for a frequency 
corner ratio k or 4. The present system uses a k = 6. This corresponds to a 
£ % 0.7, a desirable number for a simple second order system. The presence 
of a dipole suggests a £ ~ 0.5, which ironically happens to be the optimum value 
for least integral squared error in the unconstrained simple second-order sys- 
tem. The increasing of effective wind stiffness constant as a function of increased 
RMS wind velocity was found to be desirable. Throughout the work, it was felt 
desirable to maintain the maximum Phase Margin constraint. 
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RETRODIRECTIVE ANTENNAS AND THE MULTIPATH PROBLEMS 
ASSOCIATED WITH THEIR USE ON THE 
PROPOSED ORBITING DATA RELAY NETWORK 

by 

P. M. Seal* 


INTRODUCTION 

This report is based on a study of the characteristics of an Orbiting Data 
Relay Network, ODRN, designed to provide two-way communication between 
ground stations and target vehicles orbiting the earth. Essential background 
material is contained in the Final Report submitted by RCA to NASA on March 
22, 1967, entitled ’’Orbiting Data Relay Network.’’ 

The two-way communication is to be provided by including a synchronous 
satellite as a link between the ground station and the target vehicle. A large 
part of the RCA Report is devoted to a study of the types of antennas which might 
be used in such a synchronous satellite, the final recommendation being that a 
self -focusing or retrodirective antenna array be used. 

Considerable time was spent in studying the literature available on retro- 
directive arrays and the first part of the report deals with a bibliography of this 
literature along with a brief description of the highlights of each reference. 

One problem in connection with the use of retrodirective arrays is the fact 
that the array will receive two signals from the target vehicle, one a direct signal 
and the other a signal reflected from the earth's surface. This problem, known 
as the multipath problem, is particularly serious when the target vehicle antenna 
is omni-directional, and the RCA report discusses this particular problem in 
considerable detail. A large part of the summer was spent in studying various 
aspects of this problem, the aim being to determine to what extent the beam 
steering of the retrodirective array would be affected by the reflected wave. 

The second part of this report deals with the results of this study and the con- 
clusions reached. 


* Norwich University, Northfield, Vt. 



BIBLIOGRAPHY ON RETRODIRECTIVE ANTENNA ARRAYS 


Books ; 

Microwave Scanning Antennas , edited by R. C. Hansen, Academic Press, 
Volume H, Array Theory and Practice, 1966. 

Chapter 1, The Theory of Antenna Arrays, by R. S. Elliott. 

This chapter was studied mainly as a review of antenna arrays in general, 
and of scanning arrays in particular. Of special interest was the section on the 
scanned array, which shows that by varying the relative phase angle between 
currents fed to the various elements in the array, a beam may be produced 
which scans over a wide area in space, the antenna array itself being stationary. 

Volume HI, Array Systems, 1966 

Chapter 5, Self-Phased Arrays, by Donald L. Margerum. 


This chapter covers the subject of self-phased arrays very extensively, in- 
cluding a detailed discussion of the principle of operation of such arrays, and a 
description of the various types of these arrays which have been investigated. 
There is also an extensive set of references at the end of the chapter. 



Periodicals : 

The references listed below are all articles which appeared in the March, 
1964 IEEE Transactions on Antennas and Propagation. This was a special issue 
on active and adaptive antennas. 

Self-Phasing Array Antennas, by Skolnik and King. 

This article describes the principle of operation of the self-phasing antenna, 
including the Van Atta Array, and discusses how noise and signals coming from 
more than one source affect the operation of the array. It also gives a number 
of applications. 

Self-Focusing Array Research Model, by Sichelstiel, Waters, and Mid. 

Describes the results obtained on a research model containing 25 subaper- 
tures at first, then extended to 64 subapertures. Results showed that good track- 
ing was obtained in the presence of noise, even when the signal-to-noise ratio 
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was only 2 db. Good results were also obtained when signals from two sources 
separated by a few degrees were picked up by the array, even when the power 
in one signal was 12 db below the power in the other. 

Electronically Adaptive Antenna Systems, by Ghose. 


A feature of this article is that it analyzes in some detail the effect of noise 
on a self -focusing antenna on a statistical basis. It also emphasizes the im- 
portance of having three mutually orthogonal antenna arrays on space vehicles 
whose attitudes can change in an unpredictable manner, thus making the antenna 
system insensitive to the direction of arrival of the signal. 


Retrodirective Array Using the Heterodyne Technique, by Pon. 


Describes a 4-element array using a single mixer and local oscillator and 
shows its advantages over the Van Atta Array. 


An Active Retrodirective Array for Satellite Communications, by Andre and 

Leonard. ” — 

One feature of this article is that it discusses the effects of phase errors 
introduced by the amplifiers, mixers and transmission lines interconnecting 
the amplifiers. A phase error of up to ±30 degrees produces a ’’worst case" 
loss of 1.2 db in the beam output. No mention was made, hbwever, of the possible 
effect on the direction of the beam that might be caused by a phase error, es- 
pecially if the error were in the same direction for all elements. 

Also of interest was the statement that it was possible to "fabricate a circu- 
lator coupled tunnel-diode amplifier weighing approximately one ounce at S-band." 
It was also claimed that less than one watt of primary power was required for a 
200-element array. 


Spherical Retrodirective Array, by Rutz-Phillip. 

This article emphasizes the fact that a spherical array will accept a wave 
from any direction in space, amplify it and return it in the direction from which 
it came, the amplitude of the re radiated signal being independent of the direction. 
It is therefore ideal for use on satellites . 
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This array makes use of a tunnel-diode frequency converter, and one state- 
ment of interest is that "one diode functions simultaneously as amplifier, fre- 
quency translator, modulator and phase inverter." No circulator is required to 
separate the incident wave from the reradiated wave. 

The elements of the array do not need to be equally spaced, but the spacing 
between elements should be below one wavelength. Each element of the array 
should be circularly polarized, but the orientation of the axes can be arbitrary. 


A Satellite Data Transmission Antenna , by Belfi, Rothenberg, Schwartzman, 
Tilley and Wills. 

This article discusses the advantages of using a retrodirective antenna array 
in a data-collecting satellite such as TIROS IV. Upon command from a ground 
station, such an antenna would transmit its information back to the ground station 
through a beam automatically pointed in that direction. The conclusion is that 
"the hybrid-phasing-matrix antenna is the optimum system for this satellite con- 
figuration in terms of performance, weight, power drain and reliability." Several 
advantages are listed, one being that "It will not respond to jamming or an incor- 
rectly coded signal." 

A Phase-Locked Receiving Array for High-Frequency Communications Use , 
by Svoboda. 

One conclusion from this article is that "a phase-locked array can, under 
most circumstances, give interference rejection comparable to that of a conven- 
tionally phased array." However, a disadvantage is that its operation is affected 
by interference and noise in the phase detector. Of possible interest in connec- 
tion with the multipath problem is Appendix n, entitled "Effect of Interfering 
Coherent Monochromatic Signal on Performance of Phase-Locked Array." 


Satellite Communications Relay System Using a Retrodirective Space An- 
tenna , by Gruenberg and Johnson. 

The array system described is considered to be quasi-passive in the sense 
that very little power is required for its operation. The total power requirement 
for a 10,000-element array is claimed to be only one-half of a milliwatt, and 
"could be derived entirely from the ground." Several modulation techniques are 
described and a fairly extensive set of references is included. 
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THE MULTIPATH PROBLEM 


As mentioned in the introduction, a large part of the RCA Final Report on 
ODRN is devoted to multipath effects. These effects are especially serious in 
the emergency mode, where the antenna on the target vehicle is omnidirectional. 
Under such conditions the indirect signal reaching the spacecraft after being re- 
flected from the earth's surface can be quite large compared with the direct 
signal. In the RCA report, several pages of Section 4.0, particularly pages 4-38 
to 4-46, and all of Appendix D, are devoted to a discussion of this problem. 

One of the initial goals in this summer's study was to determine the effects 
of multipath propagation on the beam steering of the retrodirective array. While 
this goal was not completely realized, by any means, a few things were ac- 
complished, and the remainder of this report is devoted to these accomplishments 


Maximum Angle Between Direct and Reflected Waves . 


One of the properties of the retrodirective antenna array is that when it re- 
ceives signals from two different directions, it will retransmit each signal in the 
direction from which it was received. In the case of multipath propagation, one 
of these signals is the direct signal and the other, the reflected signal. 

All derivations contained in the RCA report assume that the angle between 
the paths taken by the direct and the reflected waves is negligibly small. The 
RCA report also shows that, in general, the reflected wave is scattered and 
arrives at the spacecraft from a wide area of the earth's surface. However, it 
also shows that the energy in the reflected wave tends to be concentrated in the 
components arriving from the specular point, that is, the point at which the angle 
of incidence is equal to the angle of reflection. 

On this basis, a study was made to determine the maximum angle between 
the path taken by the direct wave and the path taken by a reflected wave arriving 
from the specular point. Figure 1 shows the multipath geometry involved, the 
angle being referred to as A 6 . 

Equations were derived from this figure, as shown on next page, and A 6 was 
calculated as a function of the angle of incidence 6 , for the case of a target ve- 
hicle assumed to be 400 miles above the earth. 

Results are shown in the lower graph of Figure 2. For this height, the 
maximum value of A0 is about 1.6 degrees and occurs when 6 = 60 degrees, or 
when the target aspect angle © = 70 degrees. 
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T - target 
S - spacecraft 
0 = center of earth 
P = point of specular reflection 
h = height of target above the earth 
d = height of spacecraft above the earth 
r = direct distance from target to spacecraft 
r ^ = distance from target to point of specular reflection 
r 2 - distance from spacecraft to point of specular reflection 
or rj + *2 - reflected distance from target to spacecraft 
a = radius of earth 

0 = angle of incidence = angle of reflection 

(H)= target aspect angle, or angle between target's local vertical ond the line T S 

Angles 0 i/CTj $ and $2 are os *hown on Figure 
S.9 ~ angle between direct distance r and reflected distance r ^ 


Figure 1— Multipath Geometry. 
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Curves of A 0 andAf vs 0 


where A 0 = angle between direct path from 

Target to Spacecraft and reflected path, 
A f = Doppler Shift, or # Fading Rate 

0 = angle of incidence ■= angle of reflection 

Height of Target, h, = 400 miles 
Height of Spacecraft, d, - 20,000 miles 
Frequency of Transmitted Signal = 2000 MHz 
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A0 
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Target assumed to be traveling at 5 miles/ second parallel to the earth 
in the great circle plane which includes Target and Spacecraft. 

Note: Both Curves go to Zero at 0 = 0°. and 0= 90° 
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Figure 2-Curves for Obtaining Maximum A 6 and Maximum Fading Rate. 



Although calculations were not carried out for other target heights, it is 
probable that the maximum value of A# \ will be roughly proportional to target 
height. 

Equations used for calculating the angle Ad in terms of 6 were obtained 
from the trigonometrical relations of Figure 1 and are listed below: 

From triangle OTP: 


Sin n Sin 6 , 

a + r 


= 6 - Ctj' 


= (a + h) 


Sin 
Sin 8 


From triangle OSP: 


Sin 


ou 


a + d 


Sin 8, 


= 6 ~ S • 


r 


2 


= (a + d) 


Sin cfi 2 
Sin 8 


From triangle OTS: 

r 2 = (a + h) 2 + (a + d) 2 - 2(a + h) (a + d) Cos ( <t> x + <£ 2 )- 


From triangle TPS: 



also 


© = 20 + A0 - a,j = 6 + + A# . 


( 9 ) 


Differential Doppler Frequency or Fading Rate, Af . 

The fading rate is defined on page 4-43 of the RCA report as the difference 
between the doppler shift in the direct signal and the doppler shift in the re- 
flected signal. The report states that "the differential doppler is typically about 
1 kHz for orbit heights up to 400 miles." No equations were shown, however, 
and no explanation was given to show how this figure was determined. 

The equations derived for obtaining A0 can be easily extended to obtain 
equations for calculating A f , and these equations are included below. 

From these equations, calculations were made to determine Af as a function 
of 0 for the case of a target 400 miles above the earth traveling at a velocity 
of 5 miles per second in the plane made by the target, the spacecraft and the 
center of the earth. The frequency of the signal transmitted from the target was 
assumed to be 2000 MHz. 

Results of these calculations are shown in the upper graph of Figure 2. For 
the example chosen, the maximum value of Af is about 8.2 kHz and occurs when 
0 is about 55 degrees, or © is about 63 degrees. Note that when the target is 
directly below the spacecraft, or 0 - 0, Af = 0 since the target’s component 
of velocity in the direction of the spacecraft is zero. When 0-90 degrees, the 
line between the target and the spacecraft is tangent to the earth, making the 
doppler shift in the direct signal just equal to the doppler shift in the reflected 
signal. Hence Af = 0 at this point also. 

Equations used for calculating A f in terms of 0 were obtained by making 
use of Figure 1 and the equations for calculating A0 • These are summarized 
below: 


dr d(>i + r 2 ) 
dt d t 


(10) 


Af xf, 

c 
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where f = frequency of the transmitted signal 

c = velocity of propagation of signal in free space 
= velocity of light. 

Assuming that the target T in Figure 1 is traveling at a velocity v parallel 
to the earth and in the plane of the figure, 


dr 

— — = v Sin © , 
dt 


( 11 ) 


where 


and 


d( r i + r 2 ) ^dr x dr 2 y d0 d<f> 

dt \d0 + dd ) dcp X dt 


d 4> v 

dt a + h 


d$ _ 1 

d<p d<£j d0 2 

dd~ + dd~ 


( 12 ) 


(13) 


(14) 


Using equations 1 and 2, and differentiating with respect to d> 


d<fi 1 a Cos 0 

= 1 

dd a + h Cos tXj 

Similarly, using equations 4 and 5, 


dcp 0 a 

— I = 1 

dd a + d 


Cos d 
Cos a 2 


By differentiating equations 3 and 6 with respect to Q , 

dr x (a + h) 
d T = Sin e 


d0j 

Cos 4>. x Sin cp. Cot. d 

1 dd 1 


(15) 


(16) 


(17) 
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(18) 


dr. 


dd 


(a + d) 
Sin 6 


dcp 2 

Cos 0, x - — - Sin 0„ Cot. 6 
2 dd 2 


Fading Bandwidth and Coherent Bandwidth 

The meaning of the term "fading bandwidth" and a discussion of how it 
differs from "fading rate" is discussed on page 4-44 of the RCA report, and also 
on page D-13 of Appendix D of the same report. Similarly "coherent bandwidth" 
is discussed on page 4-44 and at greater length on page D-20 of Appendix D. 
Further discussion of the relation between these two bandwidths is contained on 
page 4-45 of the report and again on page 4-61. 

The only reason for commenting further on these two bandwidths is to call 
attention to the fact that the magnitudes assigned to them in the RCA report are 
all based on "moderate sea roughness." Yet on page D-19 of Appendix D a state- 
ment is made that the term used as a measure of sea roughness "could vary from 
1/3 to 3 times" the "typical" value over different terrains. 

Furthermore, Equations D-19b and D-21 for fading bandwidth, B f , as given 
on pages D-18 and D-19 of Appendix D, shows that B f varies inversely as the 
termT /a, whereas the equation for the coherent bandwidth, B coh , given on 
page D-20 shows thatB coh varies directly as the square ofT /a, where T / a may 
be considered as a "roughness factor." (Perhaps ^smoothness factor" would be 
a better term, since it should be noted that the smoother the sea is, the greater 
is the termT /a.) 

To show how B, and B . vary with the factor T /o- , the curves of Figure 3 
were drawn. Consider as an example the case of a target 100 miles above the 
earth, and the target aspect angle 0 = 0, For "normal" roughness, B f is ap- 
proximately 25 kHz and B coh is about 15 kHz. However, under a "smoother" 
sea, where T /a is increased 3 times, B f drops to 8.3 kHz whereas B coh in- 
creases to about 140 kHz. On the other hand, for a "rougher" sea, where T /a 
is 1/3 of the normal value, B f increases to 75 kHz while B coh drops to 1.7 
kHz. No attempt will be made here to discuss the significance of these figures, 
but presumably they should be taken into account in attempting to predict the 
quality of transmission under varying conditions. 
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NOTES ON A MEMORANDUM FROM P. SEALMAN OF AIRBORNE IN- 
STRUMENT LABORATORIES TO P. BARRITT, NASA HEADQUARTERS 
JULY 1, 1967 


This Memorandum proposes a frequency-hopping technique to take advantage 
of the time delay between the direct signal and the reflected signal. This tech- 
nique is well described in the memorandum and a typical example of how it might 
be used is given. However, Equation 4 on page 5 of the memo, which is the equa- 
tion used for determining the time delay , appears to be in error, mainly due to 
the fact that the term h x in the equation is assumed to be a constant independent 
of the grazing angle, i p. This error leads to an incorrect value for the time de- 
lay when = 10 degrees, which in turn leads to an incorrect value for the number 
of frequencies required, N , in the example used. 

A correct expression for the time delay is derived below, making use of 
Figures 4 and 1, as well as some of the equations already derived from Figure 1. 

When the correct expression is used in the example given in the memoran- 
dum, it is found that for a target height of 200 nautical miles, or 230 statute 
miles, and a grazing angle \p of 10 degrees, 

h Q = 147 statute miles 

making t d = 274 microseconds, instead of the 440 microseconds given in the 
example. 

This in turn will make the required value of N =10 instead of 6 and make 
the hopping rate 4000 hops/second instead of 2400. It will also change the value 
of the bandwidth B shown at the bottom of page 8 from 2.4 MHz to 4.0 MHz, and 
the value described at the top of page 9 as "under one megacycle" to an actual 
value of 1.53 MHz. 

It should be emphasized that this correction does not necessarily invalidate 
the use of the frequency hopping technique, but does appear to make its use a 
little more difficult than the example would indicate. 


Derivation of Equation for Time Delay, t d > Between 
Direct and Reflected Signal. — — — — 

With reference to Figure 4, 


t 


d 


r l + r 2 ~ r 
C 


(19) 
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[I h = 100 miles 


h = 400.* 
k miles 


h - 100 miles 
h = 400 miles. 


.h = 400 miles 


Curves of Bf and B co [ 1 vs "Rough 1 
ness Factor 11 T fo 

B f = Fading Bandwidth 

B C oh = Coherent Bandwidth 

h # Height of Target above 
earth 

(H)= Target Aspect Angle 

Velocity of Target = 8000 meters/ 
sec = 5 miles/sec 

Frequency of Transmission = 

2300 MHz 





T - target 
O = center of earth 
P = point of specular reflection 
h = Tf = height of target above the earth 
r = distance from target to spacecraft 
rj= distance from target to point P 
r^= distance from spacecraft to point P 
r j + r 2= reflected distance from target to spacecraft 

P c = extension of r to a point below the target, such that P c = PT 
Tc is perpendicular to the tangent to the earth at point P. 

Td is drawn perpendicular to - 

h Q = Te = height of target above the tangent to the earth at P. 

0 - angle of incidence - angle of reflection 
\(J = 90 ° - 0 = grazing angle 

Angles G j and ^ are same as in Figure 1 , 
a = radius of the earth as in Figure 1 . 


o 
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where t d = time delay in seconds, 

c = velocity of light in miles/sec, if r x , r 2 , and r are given in miles. 

If the distance r to the spacecraft is very large compared to the target height 
h, r and r 2 may be assumed to be parallel lines. Making this assumption, Figure 
4 shows that 


r x + r 2 - r = ed = Tc cos 6 

= 2 h 0 cos 0 {20) 

= 2 h Q Sin i p . 

2h 0 

Hence t, = Sin 0. 

c (21) 


As shown in Figure 4, h 0 is the height of the target above the tangent to the earth 
at the specular point P . Hence h 0 is a function of0 . 

The relation between h 0 , h , and the grazing angle 0 may be found from 
equations already derived from Figure 1, with the help of Figure 4. From 
Figure 1 we had 


r 


l 


(a + h) 


Sin 0> x 
Sin 6 


= (a + h) 


Sin 0 X 
Cos 0 


( 22 ) 


But Figure 4 shows that 


h Q = r x Sin 0 . 


(23) 


Hence 


hg = (a + h) Sin 0 X Tan 0 , 


where, as already shown, 


. 0 X = 6 - Oj = 90° _ 0 - otj 


(24) 


(25) 
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and 


Sin a- .= Sin 9 

a + h 


2 

a + h 


Cos ip . 


(26) 


CONCLUSIONS 

The information contained in this report should not be considered as a solu- 
tion to the multipath problem, by any means. In fact, it does not come close to 
answering the question of the effect of the reflected wave on the steering of the 
beam put out by a retrogressive antenna array. It does, however, contain a 
fairly extensive bibliography of articles on this type of antenna array. It also 
calls attention to some factors over and above those which are contained in the 
RCA report on the Orbiting Data Relay Network, factors which should receive 
careful attention in a continued study of the multipath problem. In this report, 
it is hoped that it will be of some value. 
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TRANSFER FUNCTIONS OF THE ROSMAN I ANTENNA STRUCTURE 

by 

T. G. Toridis* and Yu Chen** 


The purpose of this project is to determine analytically the transfer functions 
associated with the structural part of the Rosman I Antenna. The antenna control 
system consists of various blocks which are interconnected to form a servo 
loop. A very simplified representation of the system is illustrated in Figure 1. 



Figure 1— Simplified Representation of the Servo Loop 


Each of the above blocks may be subdivided into a number of subsystems to 
obtain more sophisticated models representing the antenna system. If the trans- 
fer functions for each block or subsystem have been computed they may then be 
used to determine the characteristics of the control system. A detailed exami- 
nation of the antenna structure indicates that it can conveniently be subdivided 
into several substructures. This is illustrated in Figure 2, on the following page: 


*The George Washington University, Washington, D. C. 

**Kutgers, The State University, New Brunswick, New Jersey. 
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Figure 2 — Simplified Representation of the Structural Loop 


Each substructure may be replaced by one or more points or lumped masses 
in a simplified model representing the entire structure. The location of these 
points may be chosen to coincide with certain key locations of the antenna such as 
the drive system, the X and Y bearings, the reflector dish, the feedbox, etc. 

In this manner, results obtained by using the theoretical model can be compared 
with experimental values representing direct measurements obtained by means 
of recording apparatus mounted on the antenna structure. The points selected 
for this purpose are listed in Table I, below: 

Table I 

Lumped Mass Points 


Points 1 and 2: 
Points 3 and 4: 
Point 5: 

Points 6, 7, 8, 9: 
Point 10: 

Points 11 and 12: 
Point 13: 


X-Wheel Bearings 
Hydraulic Drive 
X-Wheel Counterweight 
Reflector Dish 
Y -Wheel Counterweight 
Y- Wheel Bearings 
Feedbox 


The mass of the structure is lumped on each of the above 13 points. Con- 
sidering a translational degree of freedom along each of the x, y , and z direc- 
tions of a Cartesian coordinate system a simplified model of 39 degrees of free- 
dom is obtained (the x andy directions coincide with the directions of the X-wheel 
and Y-wheel shaft respectively, and z represents the vertical direction). The 
masses associated with each degree of freedom are given in Table n (M = diagonal 
mass matrix in lb-sec 2 /in.). 
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Table II 

Numerical Values of Lumped Masses (lb-sec 2 /in) 


m 4,4 = m 5,5 =m 
m 7 , 7 ~ m 8 , 8 = m 
m i0,10 = m ll ,11 
m l3 , 13 = m 14 , 14 
m l6,16 “ m l7,17 
m i9, 19 = m 20,20 
m 22,22 = m 23,23 
m 25,25 = m 26,26 
m 28 , 28 = m 29,29 
m 31 ,31 = m 32 ,32 
m 34,34 = m 35 , 35 
m 37,37 = m 38 , 38 


3,3 = 666 


6,6 = 666 


9,9 


= 252 


~ m 12 , 12 
= m l5,15 
= m l8,18 
= m 21,21 
= m 24 , 24 
= m 27,27 
= m 30,30 
= m 33,33 
= m 36,36 
= m 39,39 


= 252 
= 1200 


= 88 
- 88 


= 88 
= 88 
= 300 
= 156 
= 156 


= 96 


In order to simulate the actual behavior of the structure by the simplified 
model the various masses would have to be interconnected by springs in a com- 
plicated manner. The determination of the proper connections is a rather difficult 
matter if a close analogy at the higher frequency range is desirable between the 
actual structure and its simplified model. Since the static properties of the 
actual structure has been previously determined by the Martin SB038 program 
it is desirable to use this available and realistic information somehow. This is 
found possible by appropriately choosing such key locations on the structure that 
coincide with existing grid points used in the SB038 program. A small program 
is written to generate from the tape which has the complete influence coefficient 
matrix of the structure (645x645) a 39x39 influence coefficient matrix for the 
simplified model used in this study. This is accomplished by applying successively 
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a unit load at each degree of freedom and letting the program compute the 
resulting deflections at all 39 points. 

The influence coefficient matrix obtained in this manner is not found to be 
perfectly symmetrical. However, in order to be able to use a simple eigenvalue 
routine available in the Goddard computer library, the influence coefficient 
matrix has been forced into a symmetrical form. In addition, because of the fact 
that the diagonal terms of the matrix do not always predominate over the off- 
diagonal terms, it has been found necessary to re-estimate some of them to 
improve the conditioning of the matrix. The resulting matrix is given below in 
the form of successive rows, row 1 through row 39, where D is used to denote 
the influence coefficient matrix. ( 

Having determined the mass and influence coefficient matrices of the sim- 
plified model, the differential equations governing its dynamical motion may be 
written in matrix form as 

DMX + X = DQ W 

where X and X represent the displacement and acceleration vectors, respectively, 
and Q denotes the force vector. Premultiplying Eq. (1) by the transpose of the 
mass matrix, M T = M, and letting A = MDM (A is a symmetric matrix) the re- 
sulting equation is 


AX + MX = MDQ 


( 2 ) 


FREE VIBRATIONS 

The natural frequencies and modes of vibration of the system may be ob- 
tained by setting the right hand side of Eq. (2) equal to zero, i.e., 

AX + MX = 0 


Let 


X = W i sin co. t ( 4 ) 

where ox and W ■* represent the j-th natural frequency and mode of vibration, 
respectively. Substituting Eq. (4) into Eq. (3) and cancelling out sin cot 

-a)2 AW j + MW j =0 (5) 
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Premultiplying Eq. (5) by (W J ) T , the transponse of W J , 

-CO 2 . (W j ) T AW j + (W j ) T MW 1 = 0 


or 


(W j ) T MW j = o)2 (W j ) T AW* (6) 

The orthogonality relationship of the modes may be expressed as 

(W j ) T MW k = 0 j / k 

Eq. (6) may be put into a more general form by letting W represent the modal 
matrix (a matrix whose j -th column represents the j -th mode) and letting fi 
represent a diagonal matrix containing the squares of the eigenvalues . Thus, 

W T MW = f2W T AW (8) 



FORCED VIBRATIONS 

The solution to Eq. (2) may be assumed to be of the form 



j 


(9) 


in which 4 >. represents a scalar function of time associated with the j — th mode 
and $ represents a column vector containing all the <£ j ’s. Substitution of Eq. (9) 
into Eq. (2) yields the following equation 

AW<D + MW$ = MDQ ( 10 ) 


Premultiplying Eq. (10) by 

W t AW$ + W t MW$ = W t MDQ (11) 


Let, 

K = W T AW 
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It can easily be shown that K is a diagonal matrix. Premultiplying Eq. (11) by 
K _1 , i.e., the inverse of K , the result is 

$ + K' 1 W T MW4> = K- 1 W T MDQ 

Recalling that A = MDM, noting that K' 1 = W" 1 M” l D~ 1 M~ 

(8), the above equation takes the form of 

$ + ft $ = K" 1 W T MDQ ( 14 > 


(13) 

1 (WY 1 , and using Eq. 


or 


$ + ft$ = W 1 M" 1 Q 


(15) 


The solution to Eq. (15) depends on the nature of the force vector Q . In what 
follows the solution to this equation will be obtained using operational methods. 


TRANSFER FUNCTIONS 

Taking the Laplace transform of both sides of Eq. (15) 

$ + ft$ = W _1 M" 1 Q (16) 


Similarly, 


X = W ® 


or 


and 


4> = W 1 X 


$ = W 1 X 


(17) 


(18) 


Substituting Eqs. (17) and (18) into Eq. (16) 

W 1 X+ftW* 1 X = W -1 M -1 Q 


(19) 
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Assuming zero initial conditions and using the properties of the Laplace trans- 
forms concerning higher derivatives, Eq. (19) may be rewritten as 

W 1 ts 2 X] + QW- 1 X = W 1 M- 1 Q (20) 

where s represents a scalar function of a complex variable. Let a diagonal 
matrix Z be defined as 


(s 2 + co 2 ) 0 0 0 

0 (s 2 + co 2 ) 0 0 

Z = • (21) 


0 0 0 • • • • (s 2 + co 2 ) 


in which, as before, o> . represents the j - th natural angular frequency of vibra- 
tion. With the use of matrix Z Eq. (20) becomes 

Z W" 1 X = W _1 M* 1 Q 

or 

X - WZ" 1 W- 1 M- 1 Q (22) 


For computer programming purposes it is desirable to avoid finding the inverse 
of the modal matrix W . Noting that 

W 1 = K" 1 W T MDM 


Eq. (22) may be rewritten as 


X = W Z" 1 K- 1 W T MDQ 


(23) 


Matrices Z and K being diagonal, their inverses can be found easily. 
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The transfer functions for the system under study may be found from Eq. (22) 
or Eq. (23) by applying a series of impulses one at a time at each of the degrees 
of freedom. Since the Laplace transform of a unit impulse function is unity, for 
each application of a unit impulse the matrix Q represents a column matrix 
containing all zero elements except a single non-zero element equal to one, at 
the appropriate position. To obtain the entire set of transfer functions it can be 
shown that Q has to simply be replaced by the identity matrix, I . Therefore, 
denoting by T the matrix containing the transfer functions 

T = WZ' 1 W' 1 M' 1 Q = WZ" 1 K" 1 W T MD (24) 


Substituting Eq. (24) into Eq. (22) 


X = TQ ( 25 ) 

Eq. (25) represents the solution to the forced vibration problem (see Eqs. (2) 
and (16)). If the Laplace transforms of the forcing functions have been found, 
the operational displacements can be computed using the above equation. The 
physical displacements can then be determined by finding the inverse transform 
of X. 


NUMERICAL RESULTS 

Based on the above analysis a computer program has been developed to 
calculate the natural frequencies and modes of vibration, and the transfer matrix 
of the system. The natural frequencies and modes of vibration obtained in this 
manner are listed on following pages . 

The above eigenvalues represent the co? , or, the squares of the angular fre- 
quencies. The corresponding linear frequencies, the f . 's, in cycles per second 
may be found from 


f . 

j 



(26) 


A simple computation shows that the above eigenvalues correspond to a frequency 
range of about 1.03 cps to 43.5 cps. This frequency range in general, shows some 
agreement with the experimentally measured frequency range. 

The transfer coefficients obtained by using the computer program are also 
shown on following pages . 
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For the above transfer coefficients, the point of application of the unit im- 
pulse is indicated by referring to the "degree of freedom number." Following 
this designation, the coefficients associated with each degree of freedom are 
given sequentially (by row). Thus, for example, the 39 numerical values listed 
after the "Degree of Freedom No. 1" are the transfer coefficients for each de- 
gree of freedom due to a unit impulse applied at degree of freedom No. 1. If 
the transfer coefficients are divided by s 2 + w 2 they will then yield the desired 
transfer functions. However, the transfer functions obtained in this manner 
represent only approximations to the actual values, since they are based only 
on the contribution due to the fundamental mode. However, the contributions of 
all 39 modes have been obtained by means of the computer program. For the 
purpose of brevity, they are not listed in this report, but may be obtained through 
the Antenna Systems branch, at Goddard Space Flight Center. 

It may be observed that, as expected, the transfer coefficients possess sym- 
metrical properties since the transfer matrix of a linear system is symmetrical. 


CONCLUSIONS 

A simplified model has been used for determining the dynamical character- 
istics of the Rosman I Antenna. The model consisting of 39 degrees of freedom 
yields natural frequencies encompassing the experimentally observed frequency 
range. The transfer functions obtained by means of the model may be used in a 
control system simulation of the antenna. More details concerning this study and 
the computer results are furnished in a paper by the same authors, reported in 
the 1967 Goddard Summer Workshop Report. 
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SPATIAL LIGHT MODULATION EMPLOYING A 
FERRO-MAGNETIC FLUID 

by 

J. P. Uldrick* and P. A. Lux** 


ABSTRACT 

This report describes some work done during the summer of 1967 
concerning the feasibility of synthesizing a fluid light modulator by using 
a ferro-magnetic fluid for transducing an electrical signal into a cor- 
responding spatial transparency distribution. Several cells using 3" x 5" 
glass plates with 0.01" to 0.10" spacings were constructed. A clear fluid 
was passed through the region between the plates at about 1 foot per sec- 
ond. In the middle of the stream a small stream of opaque fluid was in- 
serted in the clear carrier fluid to form a pattern between the plates. 
Since the flow is laminar the pattern will hold its shape as it traverses 
the cell. The storage time is adjusted by regulating the flow rate of the 
clear carrier fluid. 

A ferro-magnetic fluid was manufactured for useas the opaque fluid. 
Since this fluid is highly responsive to magnetic fields^ then it is possible 
to deflect it with an electromagnet. This gives a method for transducing 
an electric signal into an optical pattern. By properly shaping the poles 
of the electromagnet it is possible to exert a constant body force on the 
ferrofluid in some region between the poles. 

It was necessary to design and construct an attritor for manufacturing 
the ferrofluid. The ferrofluid consists of Fe 3 0 4 particles of about 100 
antroms in diameter suspended in JP-4 jet fuel. A normal ball mill is not 
capable of grinding the iron powder small enough. 


INTRODUCTION 

The object of the work reported herein was to investigate the feasibility of 
storing information in a moving fluid. Ultimately the information was to be read 
out using a laser. 


*United States Naval Academy, Annapolis, Maryland 

**Sacramento State College, Sacramento, California 


185 



A cell was constructed which consisted of two 3" X 5" parallel glass plates 
spaced close together and a fluid was passed between the plates at relatively low 
velocity. Because of the close spacing and the low velocity the fluid flow is 
laminar and has uniform velocities in planes parallel to and between the glass 
plates. One such cell is shown in Figure 1.1. 

If a clear carrier fluid is used between the plates and a pattern of dark fluid 
is immersed in the carrier fluid between the glass plates near the upstream end 
of the cell its pattern is carried downstream through the cell with very little 
distortion. The velocity with which the pattern traverses the cell can be adjusted 
by regulating the flow rate of the clear carrier fluid. In this way the storage 
time can be controlled. 

Although this investigation did not progress that far, the ultimate use for this 
cell would be in an optical data processing system. In this application a laser 
beam would pass through the cell normal to the glass plates and the pattern of 
dark and clear fluid in the cell would be detected with an optical filter. 


The fluids used in the study were JP-4 jet fuel for the clear fluid and a ferro- 
fluid for the dark fluid. The manufacture of the ferrofluid will be discussed below. 


THE FLUID FLOW PROBLEM 

The principle that will be discussed to transport the desired transparency 
across the aperture of the optical system will be the steady laminar flow between 
parallel boundaries. In order to use such a system one must be able to control 
precisely the trajectories of any given fluid particle. This requirement precludes 
any possibility of local fluctuation in the velocity or mixing of the fluid as it 
passes between the boundaries. Such will be the case if the characteristic flow 
parameter composed of the product of the mean velocity V, density p, and the 
spacing between the boundaries L, divided by the dynamic viscosity p of the 
fluid (i.e., a Reynolds’ number), has a magnitude not greater than a certain 
limiting value which is roughly of the order of 10 3 . For a fluid such as kerosene 
with a kinematic viscosity p/p = 2 x 10" 4 ft 2 /sec requires that 

V L <0.2 ( 2 ‘!) 

where V and L are defined above. 

Let us now consider the laminar flow of an incompressible fluid between 
infinite parallel plates as shown in Figure 2.1. 
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The equation of motion for an incompressible viscous fluid is 


BV BV BV BV 

X X „ r X X 


Bp 


p — -+V — — + V — i +V — -:cB + Iii + uV 2 V 
H Bt x Bx y By . z Bz H * Bx ^ 


BV BV BV BV 


BV_ BV. BV. BV. 

Bt ’ x Bx ' y By ‘ z 3z ^ z + Bz 


P -d- + V. + V — ^ + V, = p B. + |& + n v 2 V, 


( 2 . 2 ) 


where 


P 

{V , V , V } 

1 x’ y’ z 

P 

{B , B , B } 

• x* y ’ 7 . 

p 


- fluid density 

- velocity components 

- pressure 

- body force components 

- dynamic viscosity 


For the case under study 



^v x 


_BV z 

Bt 

Bt 

Bt 


If the body force is conservative, i.e., B = VB where B is the body force potential, 
then the equations given in (2.2) simplify to 
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^ OoB + p) = 0 


(2.3 a) 


J-0oB + p) = O (2.3b) 

By 


_a_ 

Bz 


B 2 V 

OB + p) = n - 

By 2 


= JB 


(2.3c) 


Equation (2.3c) can be integrated to give 


y z 


= /3 (f_ 

li \ 2 


+ c x y + c 


■) 


(2.4) 


where C x and C 2 are constants of integration which can be evaluated by imposing 
the no slip condition at the boundary y = ± h. By substituting these boundary 
conditions into Equation (2.4) yield 


V = JL (h 2 - y 2 ) 

2 2/x v y ' 


(2.5) 


which describes the familiar parabolic profile depicted in Figure 2.1. If we 
orient the plates such that the xz plane is horizontal, Equations (2.3a) and (2.3b) 
specify that the pressure is independent of x and y giving no particular useful 
information here. 

Now if one introduces a small quantity, say a cubical volume of side dimen- 
sions kh (k < 1) of a dyed fluid at the point x = y = z = 0, then at any subse- 
quent time t later the center of the dye will be located at 




( 2 . 6 ) 


190 



according to Equation (2.5). The particles on the surface of the cube which was 
originally located at z = 0, x = 0, and y - kh will be located at 


z =^1(1 -k 2 )t (2.7) 

2/x 


any time t later. Thus we see that the size of the dyed particle relative to the 
plate spacing is a measure of the smear along the z axis. In the present appli- 
cation k should be made as small as possible. In principle then one should be 
able to calculate the dye distribution and hence the transparency distribution of 
the introduced dye as a function of space and time. 

In order to use this principle for spatially modulating light a method must 
be devised to introduce the time dependent electric signal into the fluid. For this 
we have synthesized a magnetic responsive fluid similar to that employed by 
Papell and Faber of the Lewis Research Center and described by Neuringer and 
Rosensweig. The production and properties of this unique fluid will be discussed 
in Section IV below. 

Basically the ferrofluid is a colloidal suspension of very fine (sub-micron) 
sized particles of ferrite such as Fe 3 0 4 immersed in a carrier fluid such as 
JP-4 jet fuel. These fluids are highly responsive to magnetic fields. According 
to Neuringer and Rosensweig the force per unit volume induced into the fluid in 
the present of a magnetic field is given by 


F/AV = /x 0 MVH * 2,8 ' 

where is the free space permeability, M is the magnetic moment per unit 
volume of the magnetic field, and H is the magnitude of the magnetic field. For 
magnetic fluids which are perfectly soft in the magnetic sense, the magnetization 
M vanishes when the applied field H is zero while at very high applied fields the 
magnetization must approach the constant, saturation value. For low to moderate 
applied fields it is reasonable to assume an empirical relationship of the form 

M = X H ( 2 - 9 ) 


to. hold, where the susceptibility is taken to be of the form 

n-2 


X = cr (T) H 


( 2 . 10 ) 



in which a (T) depends upon temperature only. Thus cr and n are magnetic ma- 
terial properties which must be determined by experiment. 

By substituting Equations (2.9) and (2.10) into Equation (2.8) and by assuming 
an isothermal exchange of energy the force per unit volume can be written as 



F/AV = V 


l- L o °" H n 


n 


( 2 . 11 ) 


This indicates that the magnetically induced body force is conservative where the 
potential is proportional to the nth power of the magnitude of the magnetic field. 

With the above results in mond let us study the following problem. Suppose 
a particle of opaque ferrofluid synthesized with a kerosene fluid as carrier be 
introduced into a transparent kerosene fluid which is flowing in the z direction 
between parallel glass plate as shown in Figure 2.2 and as discussed above. 


♦ ttt ♦ 


tttt 

— FERREFLUID PARTICLE 


1 




a MAGNETIC FIELD REGION 


Figure 2.2 
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As the particle moves along the z -axis with constant velocity V z , it moves 
through a time dependent magnetic field which is generated by an electromagnetic 
external to the flow region. The potential for one convenient magnetic field can 
be written as 

4> = [A + B cos (cot +<//)] y (x + c) [U_ 1 (z) -U_ 1 (z - a)] (2.12) 

in which K x is a measure of the strength of the field, B < A are electrical modu- 
lation constants, C is a constant which depends on the shape of the poles, and 
U_j (z ) is the unit step function and 0 is a phase angle. See Section III below for 
vari types of field shapes . In terms of the potential 0 the magnitude of the 
magnetic field is given by 



By substituting Equations (2.12) and (2.13) into Equation (2.11) it can be shown 
that 


F X /AV = jj-q cr n {Kj 


[A + Bcos (cot +0)] n [U 


(z) - U_j (z -a)] ,{y 2 + (x +c) 2 } n/2 “ 1 (x +c)} 


F y /AV = ^ orn.iKJ [A+Bcos (cat +0)] n [U_, (z) -U_ t (z - a)] y 

F z ~ 0 


(2.14) 


If the particle entering the magnetic field is traveling in the x - y plane 
then according to Equation (2.14) it will experience a force 

F X = cr AV n Kj (A + B cos (cot +0)) n (X+C)"" 1 [U_ x (z) - U_ x (z-a)] 

F v-° 

F = 0 (2.15) 

Z 


while in the magnetic field region. 
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Now it is desirable to be able to predict the location of the particle at any 
subsequent time after it is past the region of the magnetic field. By assuming 
laminar flow of the transparent fluid around the magnetically energized fluid 
particle, the drag experienced will be directly proportional to the relative 
velocity between the ferrofluid particle and the surrounding transparent fluid, i.e., 

drag - Dx 

where D is a constant depending upon the shape and size of the ferrofluid particle 
relative to the plate spacing and the viscosity of the transparent fluid, and x its 
velocity in the x direction. Then, to a very good approximation, the equations 
describing the motion of the ferrofluid particle will be 


p { A V x = AV fx Q cr n K £ (A +B cos (cot +'/'))" (x + c) n_1 [U_ 1 (z) - U_ 1 (z - a)] - Dx 

y = 0, z = 0 (2.16) 

where p f is the density of the ferrofluid. The initial condition are at t = 0: 
x=y=z=x=y=0, and z = V z = const. 

As pointed out above the exponent n is a property of the magnetic fluid. 
Neuringer and Rosensweig reports a value of n = 4/3 for their fluid. It is rea- 
sonable to expect that the ferrofluid we synthesized at Goddard will possess mag- 
netic properties of the same order of magnitude. If this is the case, then the 
forcing function given in Equation (2.16) depends rather weakly upon x, i.e., x 1/3 . 
Furthermore, by assuming the induced displacement of the particle to be small 
while moving through the region of the magnetic field we may simplify our anal- 
ysis considerably by saying that 


Cj = (x+c) n_1 (2*17) 

With this simplification Equation (2.16a) can be written in the form 

x + ax = 0 [A + B cos (a>t + '/')] n [U_ 1 (z) - U_ 1 (z - a)] (2.18) 


where 


a 


D 

Pf AV ’ 




P Q o- n KJ c x 
_ 
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To obtain the general solution of Equation (2.18) we first expand the forcing 
function in a Fourier series as 


/3 [A+B cos (cot + <//)]" = /3 


A 0 + 2 


00 



A m cos m (co t +\p) 


(2.19) 


in which the A can be determined from the Fourier inversion formula 

m 


A 


m 


— I [A+B cos (cot + i//)] n cos m (cot + i//) d (cot +ip) 

7T I 

Jo 


By employing Laplace transform techniques, Equation (2.16) may be integrated 
to yield 


x = A 0 | — + — e- at - Z 

V a a 2 a 2 / 


00 

L 


A Gosim b 

m f 

m= l ( mw ) 2 + a 2 


, a 

cos mcot sin m cot 

mco 


-f Z ^ A n m a; sin mi// 

m=l 


m 2 co 2 (a 2 •( m 2 co 2 ) 


.(a cos m^t - majsin mcot) 


1 1 

+ 


-at 


am 2 co 2 a (a 2 + m 2 co 2 ) 


(2.20) 


for the transverse displacement of the particle, which is only applicable while 
the particle is in the region of the magnetic field. 

Thus, as the particle leaves the magnetic field, z = a , its displacement 
will be 
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x =a 0 (jl a,-'’' 1 
■ 0 l“V, a 2 


j_\ 2 y’ A . c ° sm ^ 

» 2 / ^st ( m “) 2 + “ 2 


a a . a 

cos mco — - sinmo) — 

V mco V 


uu 

+ 2 y A n mo> sin riup 

m=si 


m 2 co 2 (a 2 + mo) 2 ) 


a . a 

a cos mw— - ma> sin mco — 
2\ \ V. V.. 


1 1 

+ 


-a a /V 


a m 2 a> 2 a (a 2 + m 2 a> 2 ) 


(2.21) 


At this instant its velocity will be 


y-i A m cosm^ p 

4l* m2 " 2 + a2 L 


a a 


A o 1 - a a/v X * r . a a 

-X n— - — e +2 / |mo> sin mco—- + a cos mco 


00 

-2 y A m ma> sin m\p 

m= 1 


1 / . H 9 2 3 

- araw sin m a) — + rrr co^ cos mco — 

|_m 2 a> 2 (a 2 +m 2 a> 2 ) \ z> 


1 -a a/V, 

+ e 

a 2 + nv 2 co 2 


( 2 . 22 ) 


To determine the transverse location of the particle at any subsequent time 
after leaving the magnetic field, we have for the motion equation 


x = - ax 


(2.23) 


This can be integrated simply to give 

-a(t— a/V ') 


x = x e 

a 


(2.24) 


and 


x =x a - 


a 




(2.25) 


196 



The z displacement of this particle will be 

z — V t (2.26) 

Z 

Therefore, a particle location at z = V z t will have a lateral displacement 


x 


= X 


a 



(z/V z 


a/V, 


>] 


(2.27) 


In order to calculate the shape formed by a continuous line of streaming 
particles, it is necessary to relate the phase angle 0 of the magnetic field with 
the particular particle under consideration. Assuming the phase angle of the 
magnetic field was 0 O when the particle which is located at z = a entered the 
magnetic field at z = 0, then the phase angle for the particle located at any z > a 
would be 


0 = 0 Q + a) (2.28) 


By substituting Equation (2.28) into Equations (2.21) and (2.22) and then 
substituting the results into Equation (2.27), we can derive the shape of the re- 
sulting pattern of streaming particles at any instant. Due to lack of time, we 
were unable to pursue the analytic solution further. 

In order to see if this type modulator was feasible, several cells were con- 
structed and several different modulation configurations were tried. Descriptions 
of some of these will be found in the next section of this report. 


DEFLECTION SYSTEMS 

In this section we discuss several possible methods for deflecting the stream 
of ferrofluid which is being transported in the nonmagnetic clear carrier fluid. 
Each of these methods are discussed on the premise that the ferrofluid can be 
displaced in the carrier fluid by inducing an electromagnetic body force in the 
ferrofluid. 

The basic equation for governing the design of possible deflection systems 
is that given by Equation (2.11). Two of the systems studied are analyzed on the 
next page. 
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Single Wire Normal to the Fluid Flow 


As demonstrated by Rosensweig, a ferrofluid with a high magnetization M 
can have impressed upon it a measurable force by passing a current through a 
wire located near the fluid. Thus , let us consider the force induced in a small 
volume A V of ferrofluid which is located near a wire carrying a current I as 
shown in Figure 3.1. 


WIRE’ 


I'M: 


FERROFLUID PARTICLE 


Figure 3.1 


It is well known that the magnetic field in the region exterior to a current 
carrying straight wire whose center is located along the y axis is 

O 

3 = -L ? (3-i) 

277T 


where I is the current, r is the radial distance from the center of the wire to 
the point in question, and r is a unit vector normal to the radius vector r and 
the unit vector X along the y axis. By substituting Equation (3.1) into Equation 
(2.11), the force acting on a unit volume located at (x, o , z) is 




3r 3 r 

3 x 3y 



(3.2) 
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But, in the present case r = / x 2 + y 2 . Thus, equation (3.2) reduces to 


F/AV = fi 0 a 


(— 

Vi 

2x 

”7* 2z 7* 

i + k 

\ 2rrr 

/ r 

r 3/2 

r 3/2 


Assuming a value of n - 4/3 this last expression can be written as 


(3.3) 



/_I_V 4/3 xl + zk 

\ 27T ) (x 2 + y 2 ) 17/12 


(3.4) 


This equation reveals some interesting results concerning the force on the 
small particle as it passes the wire. For a small volume of ferrofluid, the 
maximum force exerted will occur when the particle is closest to the wire. For 
example, if the particle is in contact with the wire as it crosses the x axis, the 
force will be 


F = 


2 a A V 



R~ 11/6 


(3.5) 


where R is the radius of the wire. If AV is very small, either the wire must be 
very small or the current through the wire must be very high to expect any 
measurable force. However, the maximum current is limited by the burnout of 
the wire. The burnout current can be approximated by I max = KR 4 . Thus we 
see that 


2 cr A V 


KR 1673 

(2tt) 4/3 


eR 3 ' 5 AV 

r11/6 


(3.6) 


where e is a bounded constant. This last equation indicates that this system is 
impractical in the present application. 
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Deflection System Using a Shaped Pole Electromagnet 


There are several ways to shape the poles of an electromagnet for use in 
deflecting the stream of ferrofluid. Since we are interested in a system where 
the pole faces are exterior to the fluid region, we use the following analysis. To 
find the shape of the poles it will be assumed that the permeability of the core 
material is of sufficient quality such that the pole faces form magnetic equipo- 
tential surface. Since /z of iron is greater than 1000 , this is a reasonable as- 
sumption. The problem now at hand is to determine the shape of the equipotential 
surfaces to give the desired field along the x axis. 

Let the fluid be flowing along the z direction. It is desired to induce into 
the ferrofluid a force along the x axis with a magnitude which is approximately 
independent of x. Thus, according to Equation ( 2 . 11 ) we set 

F/AV = ^<rH”-‘(l*i + |JiT + |5.E) = - K? < 3 ' 7 > 

\ ox ay a z j 


where K is weakly dependent upon x. Also, the H vector field must satisfy 

Curl H = 0, div H = 0. ( 3 - 8 ) 

In order to solve Equation (3.7) along with (3.8), we assume that 17 is close 
enough to unity that /z 0 a H n_1 = constant. Also, in order that the pole faces be 
exterior to the flow we require that the H field be normal to the x-axis, i.e., 

H = H = 0 here. Thus we write 

X z 


BH _ BHy 
Bx Bx 

where 


(3.9) 


“1 7 

Now according to Equation (3.8) 

H = -V0 and V 2 <£ = 0 


(3.10) 


where <j> is the magnetic potential. 
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I 


Therefore, along the x axis 



(3.11) 


One solution that will satisfy both Equations (3.10) and (3.11) is 


4 > - Kj y (x + C 2 ) ( 3 - 12 ) 

Now to find the shape of the pole faces , it is only necessary to set 4 > equal 
to a constant and to choose a suitable value for C 2 . For example, the shape of 
the pole faces used in the experimental part of this program it was assumed that 

4 >/ Kj = ± 0.06 

c 2 = - 0.6 (3.13) 


Substituting these values into Equation (3.12) the shape of the poles becomes 


± 0.06 
(x - 0.6) 


(3.14) 


The pole shapes as given by Equation (3.14) are shown in Figure 3.2. 


DISCUSSION OF PROPERTIES AND PRODUCTION OF FERRO- 
MAGNETIC FLUIDS 

Ferromagnetic fluids are colloidal suspensions of very small (sub-micron) 
size ferromagnetic particles immersed in an organic carrier fluid such as kero- 
sene. A dispersing agent is added to prevent flocculation of the iron particles. 
These ferrofluids are not commercially available. To synthesize a ferrofluid, 
a mixture of magnetic iron powder, a surfactant, and the carrier fluid is prepared 
and ground for several days in a ball mill. 

Papell and Faber of the Lewis Research Center synthesized a magnetic 
colloid which they employed to simulate zero- and reduced gravity to study the 
effect of gravity on pool-boiling systems. 
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Neuringer and Rosensweig discuss the fluid dynamics and thermodynamics 
of magnetic responsive fluids. They treat the fluid as a continuum and develop 
what they refer to as the Ferrohydrodynamics equations of motion. 

These ferrofluids are quite different from the magnetic clutch materials 
used in the 1940’s. The magnetic fluids discussed here are truly colloidal in 
that the sufficiently fine particles can be suspended indefinitely in a liquid even 
though the density of the magnetic particles is greatly different from that of the 
liquid. A qualitative explanation of the mechanics of interaction of the particles 
will be helpful in understanding its manufacture and its response to magnetic fields. 
If we consider a distribution of small spheres in a liquid we can envision that the 
thermal agitation of the molecules striking the individual spheres will prevent 
them from flocculating and settling out. Now should these spheres be magnetized 
the magnetic energy of attraction is proportional to the square of the radius of 
the particle if they touch and inversely proportional to the cube of the radius . 
Therefore, if the radii of the spheres are made small enough, it is conceivable 
that the random impact of the molecules would be sufficient to prevent floccula- 
tion of the magnetized particles . For normal temperatures and for the magnetic 
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materials available, it can be shown that particles of 25 - 100 Anstroms in diam- 
eter would be stable under these conditions. However, when two particles of 
such small diameter are brought into contact there will be an attraction present 
due to the Van der Waal forces. To counteract this latter type force, we can 
imagine coating each sphere with an elastic material whose force originates 
from the deformation of the surface. Therefore, it can be seen that the secret 
for synthesizing a ferrofluid is to grind the ferromagnetic particles to a very 
fine size and at the same time use a surfactant to coat the individual particles. 

2. Heptane Ferrofluid in Water Carrier. At our request, Stephen Papell of the 
Lewis Research Center, sent us a very small quantity of the ferrofluid they were 
using in some of their work. This fluid was synthesized with a Heptane carrier 
base. We tried to use this ferrofluid in place of the dye in the water cell. This 
proved to be fruitless because the heptane ferrofluid immediately balled up and 
diffused against the glass plates. 

3. JP-4 Jet Fuel and Heptane. By using clear JP-4 jet fuel for the main flow 
and heptane ferrofluid for the dye, a very fine line could be maintained. This 
line could be deflected quite readily with a permanent magnetic. 

4. JP-4 and JP-4 Ferrofluid. After some ferrofluid of our own manufacture was 
obtained, several types of deflection methods were tried. These methods are 
described below with comments about their operation. Again, the order is 
chronological. 

a. Single wire transverse to the flow . A small hole was cut through both 
glass plates and a wire was passed through the holes. This wire was 
then transverse to the fluid flow and was small enough so that it did not 
have much effect on the flow. A small stream of ferrofluid was inserted 
close to the wire. A current was then passed through the wire and the 
effect on the stream noted. 

The effect on the stream was zero. The measurement techniques 
were not extremely accurate, and this method was abandoned for the 
present. 

It is possible that this method may still work. Referring to Rosen- 
swig r s paper, it is noted that he was able to exert considerable force on 
a ferrofluid using a current through a wire. For example, he caused a 
ferrofluid to climb up a wire with a current in it. 

This method has an inherent non-linearity which makes it potentially 
difficult to control. For a theoretical discussion see Section III. 
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b. A single 0.006" thick lamination . A transformer lamination 0.006” thick 
and 1/2” wide and 1-1/2” long was inserted through the side of the cell 
parallel to the glass plates. Part of this lamination was left sticking 
outside the fluid and around this part a 100 turn coil was wound. A ferro- 
fluid stream was passed close (approximately 0.1”) to the end of the 
lamination which was in the fluid. 

Since we were unable to obtain a sufficient quantity of existing magnetic 
fluid for our purposes , it was necessary to produce a quantity of fluid here at 
Goddard. 

First, an attempt was made to grind the mixture using a standard ball mill. 
After two weeks of continuous grinding with a ball mill located in the Material 
Research and Development Branch, no appreciable colloid had formed. After 
a telephone conversation with Stephen Papell of the Lewis Research Center, we 
decided to design and construct a fast -grinding machine or attritor. The details 
of the attritor are shown in Figure 4.1. 
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The attritor was mounted on a drill press (as shown in Figures 4.2(a) and 
4.3) which can be driven at four different speeds. To make the ferrofluid, the 
attritor was charged with 450 cc of JP-4 jet fuel, 25 cc or oleic acid, and 60 
grams of Fe 3 0 4 powder along with 5000 3/16" diameter chrome alloy steel ball 
bearings. The attritor was run at approximately 490 RPM continuously for one 
week. Cooling water was circulated through the jacket for temperature control 
of the mixture. To determine if any colloid was being formed, we extracted 
a small sample of fluid from the attritor and centrifuged it. If after one day 
of grinding the sample appeared opaque and responsive to a small permanent 
magnet, we were confident that a ferrofluid was forming and the run was con- 
tinued. The response of a sample to a permanent magnetic is shown in Figure 4.4. 
This particular sample was extracted at the end of three days’ grinding time. 


Experimentation 


In order to ascertain whether this method of modulation has merit, we con- 
structed several cells and experimented with different techniques for converting 
the input signal into some type of corresponding transparency pattern. These 
experiments are discussed in chronological order. 

1. Water Cell . To demonstrate that the laminar flow could possibly be used for 
signal storage, a cell was constructed using two 3" x 5" photographic glass 
plates with a spacing of 0.012". Strips of adhesive tape was used for spacing the 
plates and two 1/4" rubber hose were split and mounted on the ends for channeling 
the fluid through the plates. The system was sealed with epoxy. Water was 
syphoned through the assembly and a food coloring dye stream was inserted into 
tiie upstream portion between the plates with a hypodermic needle. The system 
is a type of Hele-Shaw apparatus. It was demonstrated that a very fine steady 
line could be maintained. A photograph of this cell is shown in Figure 5.1. 

When a current was applied to the coil the fluid was deflected. The problem 
was that the fluid tended to stick to the lamination. This no doubt is due to the 
large magnetic gradient at the surface of the lamination. 

c. Pair of laminations outside the fluid . Two laminations of the size of the 
previous case were placed normal to the glass plates, one on either side 
of the cell, as shown in Figure 5.3. Coils were wound on them. They 
were energized to give first two north poles against the glass and then a 
north and a south against the glass. 

The ferrofluid stream tended to stick to the inside of the glass plates 
and it was not possible to deflect the stream in a useful way. 
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d. A tape wound core having a cross section of 1/4" X l/4" was cut to form 
pole shapes as derived in method a. above. The gap was about 0.2” on one 
edge and about 0.4" on the other. The shape was approximately that 
shown in Figure 3.3 and mounted as shown in Figure 5.4. 

This electromagnet was used in conjunction with a cell with 0.1" spacing 
between the plates. A 300 turn coil was wound on the torried and a d.c. 
current of 0-1 ampere was passed through the coil. 

The stream was deflected, but due to the mechanical instability of the 
system it was difficult to make any accurate measurements. For example, 
slight heating of the glass plates by the magnet caused large deflections. 
Another problem with this type of deflection is the large area over which 
the magnetic field acts . In some cases, it appeared that the magnetic 
field might cause a longitudinal instability of the ferrofluid stream. 

To continue in this direction of investigation, it would be necessary 
to construct a second generation of fluid cells where the various param- 
eters are under better control. 


CONCLUSION 

It has been demonstrated that the laminar flow of a fluid between parallel 
plates can be used to store information in the form of an opaque fluid immersed 
in a moving transparent carrier fluid. This is possible because in laminar flow 
a signal of dark fluid inserted in the clear carrier fluid holds its shape as it 
flows through the space between the parallel plates. The plates are made of glass 
and a laser light beam could be passed through the storage cell to read out the 
stored information. 

Several cells of different configurations were constructed to determine the 
problems associated with immersing information into the fluid stream between 
the plates. The method which appeared most promising was to use a small stream 
of dark ferrofluid in a clear carrier fluid. This small stream can be deflected 
by an electromagnetic. However, due to the fact that this method involves the 
mechanical movement of the dark fluid relative to the clear carrier, and also 
because of difficulties in containing a desirable magnetic field within a small 
region, the frequency response is poor. Due to these difficulties, it is recom- 
mended that other techniques be used in generating the pattern. 
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Other methods that warrant consideration are the following: 

1. Use an electromagnetic ferrofluid transducer to control the flow rate 

of the opaque fluid before immersing it in the clear carrier fluid. For example, 
wind a coil around the hypodermic needle. 

2. Design a surface wave generator employing the ferrofluid and read out 
the stored information out by optical techniques. The advantages of this method 
lie in the fact that a simulated gravity field can be generated with a permanent 
magnetic field thus allowing control of the wave propagation speed and at the 
same time it appears quite easy to generate the desired surface wave pattern 
with another independent electromagnet. 

3. Use a cell as described herein but generate the transparency pattern with 
techniques that will not disturb the main flow of the clear carrier fluid. For 
example, generate the dye pattern by a flash Photolysis process as described 

by Popovich and Hummel. 
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AUTOMATIC PROCESSING OF IONO GRAMS 

by 

William A. Barrett ' 


The project that I undertook this summer as an ASEE-NASA Faculty Fellow 
had as its purpose the automatic processing of data from the topside sounder 
experiment aboard the Defense Research Telecommunications Establishment’s 
Alouette satellite series. This data, commonly referred to as an ’Ionogram’, is 
used to study the ion density in the upper layers of the atmosphere. This is ac- 
complished by transmitting a swept-frequency signal and recording the time re- 
quired for the echo to return and also the magnitude of the echo. 

In my work I was involved with Dr. F. P. Kuhl in the formation of the logic 
structure that would serve to automatically classify the signal. This phase of the 
project consisted of modifyingthe procedure used in a digital computer character 
recognition program. The modification was necessary to overcome the noisy 
background of the ionogram. We accomplished this phase by using the main 
characteristics of the signal as a key. 

While Dr. Kuhl proceeded to write the precise logic, I turned to enhancing 
the analog signal. My method called for the signal to be repeatedly integrated 
and squelched. 1 hoped that this would decrease the noise by essentially serving 
as a sharp cutoff filter. After assembling the necessary equipment and design- 
ing and building a squelch circuit, I found that this method gave a signal-to-noise 
voltage gain of approximately two. 

This work is published in X-565-67-380. 
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EXPERIMENTAL AND THEORETICAL STUDIES 
OF THE CARBON DIOXIDE LASER 

by 

John R. Bolte t 


The Quantum Optics Section at the Goddard Space Flight Center is con- 
cerned with the problem of developing a deep space communications capability 
using a high powered carbon dioxide laser operating at 10.6 microns. Part of 
this work involves a careful study of the characteristics of the laser and is under 
the direction of Mr. Nelson McAvoy. Specific studies during the summer have 
included: 

1. Lifetime Studies of Sealed Off Laser Tubes . Until very recently C0 2 
lasers had been operated only in gas flow systems because the lifetime of sealed 
off tubes was very short. Tube lifetimes of 2000 hours now appear to be possible 
using CO 2 and helium mixtures and careful tube cleaning techniques. 

2. Improvement of Techniques for RF Excitation of Plasma Tubes. The 
problem is to find efficient means of coupling rf power into the laser tube. New 
rf electrode designs are now in use which confine the electric field to the plasma 
region and additional studies are to be made in which the electrodes are placed 
inside the glass cooling jacket. 

3. Experimental Work with Detectors for 10.6 Micron Laser Radiation. 
Thermal detectors in the form of thermopiles and bolometers as well as semi- 
conductor detectors are available in the 10.6 micron infrared region. Thermal 
detectors have the disadvantage of having a slow response time while semi- 
conductor detectors have to be cooled to liquid nitrogen or liquid helium tem- 
peratures. All of the devices have good sensitivity, however. Thermistor bridge 
detectors were constructed during the summer which could easily detect infrared 
radiation from the hand held several feet from the detector. 

4. Theoretical Studies of the Resonant Modes in a Fabry-Perot Optical 
Cavity. The journal articles by Fox and Li, Boyd and Gordon, and by Kogelnik 
have given abbreviated theoretical treatments of optical resonant cavities. A 
thorough understanding of these articles is necessary so that meaningful experi- 
mental work can be done in the design of mode stable lasers. Further theoretical 
study as well as experimental work is presently underway which should lead to 
the development of highly stable lasers which will be suitable for communications 
work. 
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A PERTURBATION THEORETIC TREATMENT 
OF THE FERMI FUNCTION 

by 

B. Chern + 


During my tenure as an ASEE-NASA Fellow at the Laboratory for Theoreti- 
cal Studies of the Goddard Space Flight Center, I launched a Field-theoretic 
perturbative treatment of the Fermi function which is used in nuclear beta decay. 
A purpose of the investigation was to see whether the logarithmic divergence 
which occurs in the Fermi function when the nuclear radius goes to zero also 
occurs in the perturbative treatment. The larger overall objective is to obtain 
precise values of the vector coupling constant and to see whether the radiative 
corrections to beta decay can be made consistent with a Universal Fermi inter - 
ation. 

This work has direct application to nuclear beta decay processes in stellar 
interiors 
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BIT ERROR PROBABILITY TESTS FOR MEASURING THE 
PERFORMANCE OF THE UNIFIED S-BAND GROUND STATIONS 

by 

David A. Dalgleish ^ 


Project: 

This summer project consisted of learning more about modern communi- 
cations theory and such basic concepts as probability of bit error, noise, digital 
codes and formats, bi-phase and phase modulation, phase detection and phase- 
locked loops, and system signal-to-noise ratios. 

Specifically, the study was directed at a Bit Error Probability test that is 
part of the Station Readiness Tests of the single 30 foot, split 30 foot, and 85 foot 
unified S-band (USB) stations that are part of the Manned Space Flight Network 
found around the world. It was desired to theoretically describe the test. The 
results are partly descriptive and partly theoretical. 

Results: 


The test uses a simulator to put out a pulse code modulated (PCM) telemetry 
signal similar to that put out by the Apollo Spacecraft. This signal has a non- 
return to zero-level (NRZ-L) code format for the PCM which is bi-phased modu- 
lated on a 1024 KHZ sub-carrier. The sub-carrier is in turn used to phase 
modulate the main 2287.5 MHZ spacecraft down-link channel main carrier, which 
in this case is a simulated main carrier. 

The simulated signal is then the input to the USB system either through the 
collimation tower and antenna or directly into the parametric amplifiers with the 
antenna held at a "quiet sky" or zenith position. This input signal goes through 
the processing of the double super-heterodyne S-band receiver. The 1024 KHZ 
sub-carrier is stripped off and, in turn the PCM data is removed in the Signal 
Data Demodulator System. The noisy serial bit train is processed through a sig- 
nal conditioner and compared with the original signal in a comparator. The bit 
error count, actually a not-error count , is made for input signals having a variety 
of signal to noise (carrier to noise) levels. 
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In the test of the USB and PCM equipment, it was felt that the theoretically 
best result was 


P. = ~ P(l/D) + jP(D/l) 


where 


P(l/D) = P(D/1) P e =- 1 

2 



where S/N is a power ratio and erf 


X 



e ~ y dy 


The actual signal to noise ratio for the various parts ofthe system can then 
be determined using effective bandwidths and temperatures. With this KTB con- 
cept, a final adjustment in the S/N can be made so that a final curve of proba- 
bility of error vs input carrier to noise ratio can be made. This curve shows the 
best possible condition in the system. Somewhat arbitrarily, an acceptable mar- 
gin of +5db was selected; i.e.,the carrier-to-noise vs probability of error curve 
can be shifted +5db to give an upper acceptable bound for the test. 

Data which has been collected from many stations in the MSFN shows this 
test to be very useful in determining system performance and finding problems 
in the loop. 

Future Work: 

Future work should be spent on diminishing the 5 db acceptable margin. 
Since the theoretical curve is obtained assuming an equal probability of "ones” 
and "zeros,” some work will have to be done on simulation signals that do not 
have equal numbers of ones and zeros. Another aspect to consider is the trans- 
ition density percentage. This means that alternate 1-0 formats will give a 
maximum number of transitions (100%) and a long series of one’s followed by a 
long series of zeros (for example 64 ones then 64 zeros) will give a very small 
number of transitions (about 2%) . The effects of transition density on probability 
of error then remains to be determined. Other probable causes of error would 
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be non-gaussian noise, a quiet sky not being quiet, inaccuracies in setting in the 
input signal to noise levels for the test, and long time temperature change ef- 
fects since the overall test takes 3-4 hours. 
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of the Manned Flight Operations Division who made this summer so meaningful. 
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AN OPTICALLY INDUCED FREE CARRIER LIGHT MODULATOR 


by 

t 

Dr. Carl L. Gruber 


My work as an ASEE-NASA Summer Faculty Fellow involved theoretical and 
experimental studies of a wideband optical modulator to be used in infra-red 
communication systems. The basic concept evolved from prior optical modulator 
development by my NASA colleague, Mr. William Richards and the theory from 
my investigation of problems related to other free carrier optical modulators. 

A theoretical model for the modulator using a semiclassical approach was 
developed during the initial investigations. Using this model the optical trans- 
mission characteristics, depending upon several material properties and modu- 
lation power level, were determined. The final phase of development involved 
experimental verification of the model. Since limited time was available, com- 
plete experimental definition of the modulator was not possible. However, work 
will continue in this direction both at NASA and at South Dakota School of Mines. 

The basic modulator is described by two components. A slab of semi- 
conducting material possessing good transmission properties at the wavelength 
to be modulated and a source of readily modulated optical radiation at the band 
gap energy of the semiconductor material (a current modulated diode laser is 
ideal). The source illuminates the surface of the slab generating an inhomogene- 
ous distribution of free carriers in the semi-conductor by electron-hole pair 
production. The generated free carriers in turn alter the complex optical index 
of refraction of the material in proportion to the incident source power level. 
Either phase of amplitude modulation may be effected depending upon material 
properties and incident source power. 

This device exhibits characteristics superior to those of other modulators 
presently in use. The bandwidth should approach 1 GHz with 10% modulation 
power efficiency (three orders of magnitude better than previous devices) and 
further requires no direct electrical connections to the modulation element. 

This work has resulted in a patent disclosure and NASA Tech Note co- 
authored by Mr. Richards and myself. A paper is to be submitted for formal 
publication in Proceedings of the IEEE pending further experimental definition 
of the device. 
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AUTOMATIC PROCESSING OF IONOGRAMS 


by 

Frank P. Kuhl t 


The ionograms obtained from the ALOUETTE topside sounder are analyzed 
to show the variation with height of the electron density in the ionosphere. Iono- 
grams can be collected over all regions of the earth for any time of day since 
the orbit of the ALOUETTE satellite traverses each region periodically every 
three months. The typical appearance of an ionogram is that of a pair of inter- 
woven lines shaped like an "S" rotated 90 degrees counterclockwise. Figure 1 
shows good ionograms. 

An automatic system for machine encoding and identification of ionograms 
was proposed by the author in the document X-565-67-380, August 1966, Goddard 
Space Flight Center, as a result of last summer’s work. This summer work 
toward implementing this system will be published in a new X-document entitled 
"Automatic Processing of Ionograms." Briefly the system consists of the func- 
tional blocks shown in Figure 2. Although complete specification of the operations 
of the classification blocks (3,4,5, and 6) is not yet possible, such classification 
can now be shown to be feasible as a result of the development of FORTRAN flow 
charts for noise removal, continuity insertion and image encoding (blocks 1 and 
2). The types of classification that can be developed are suggested in X-document 
X-565-67-380. 
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Figure 2- Encoding and Identification System Block Diagram 








PCM TELEMETRY ANALYSIS 
by 

Nicholas Kyriakopoulos ^ 


The study was concerned with signal degradation in telemetry systems. The 
ultimate goal of such study would be the development of a mathematical model 
for a PCM telemetry system consisting of: antenna, preamplifier, multicoupler, 
receiver and data handling equipment. This model would be used for computer 
simulation and evaluation of telemetry systems. 

The first phase of the program involved the investigation of the behavior of 
phase-locked loops, used in the detection and bit synchronization stages, in the 
presence of noise. The system considered was a second order loop whose phase 
modulated input was corrupted by additive, Gaussian noise. An analytical solu- 
tion of the differential equation involving the probability distribution of the out- 
put signal was sought. It has, however, been shown that the output of the phase- 
lock loop under the given conditions forms a Markov process and the differential 
equation involvingthe probability distribution of such process has as of today, no 
solution. It was, therefore, decided to abandon seeking analytical solution to 
these problems and proceed along the road of computer simulation. 

The second problem investigated was the operation of a filter in the S-band 
transponder aboard a spacecraft which is part of the Range and Range Rate sys- 
tem. A question had been raised as to how closely the actual filter approached 
the claimed maximally flat time delay filter. 

A computer program to analyze the circuit was used and the computed re- 
sponse did not approach the claimed ideal one. The investigation will continue. 
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by 


Richard L. Martin^ 


Problem Statement 


The problem considered here is: 

1. Given a pulse frequency modulated signal (PFM) consisting of finite du- 
ration bursts of constant but unknown frequency distorted by clipping and 
corrupted by Gaussian noise. 

2. Develop a computer algorithm which will extract the fundamental frequency 
component of each individual burst and give an indication of the signal- 
to-noise ratio. 

Hardware Implementation 


There is an existing hardware implementation which performs the stated 
function. It consists of a comb filter bank, switching logic and a zero-crossing 
counter. 

The comb filter bank consists of 128 narrow band filter, each having a 100 
Hz bandwidth, with each center frequency separ ated by 100 Hz covering the entire 
range of possible received signals (3,6 kHz to 16.3 kHz), The received signal 
plus noise is fed into all elements of the comb filter simultaneously. As the sig- 
nal progresses in time, the outputs of the various elements built up until one 
element reaches a predetermined nominal level. At this point a low resolution 
estimate of the frequency is obtained, i.e. , the center frequency of that particular 
element. 

The output of that element is then used as the input to the zero-crossing 
counter. The total noise power contained in the signal has been reduced by the 
ratio of the filter element equivalent noise bandwidth to the input noise band- 
width. Also the frequency of the signal is known within 1%. These two results, 
a relatively clean signal and a fairly good estimate of the frequency, are neces- 
sary for the operation of the zero-crossing counter. 


^University of Maryland, College Park, Maryland 
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Methods Considered for Computer Implementation 


Several methods were considered for the digital computer frequency meas- 
urement, all of which use a zero-crossing counting scheme for the final high ac- 
curacy measurement. 

The first method, a digital simulation of the comb filters, was rejected im- 
mediately as requiring a prohibitive member of computer manipulations. Assum- 
ing a 10 ms burst with a 20 kHz bandwidth, sampling of the signal must be done 
at a minimum 40 kHz rate resulting in at least 400 samples. Processing these 
400 samples through 128 different filters was considered an excessive amount of 
work. 

A second method consists of cross-correlating the input signal with 128 dif- 
ferent square waves, the frequency of each being the center frequency of the 
filters stated previously. By using unit magnitudes for the square waves each 
cross-correlation is reduced to a summation rather than an actual multiplication 
plus summation procedure. If the input signal frequencies were from a discrete 
spectrum this would probably be the best technique to use. However, since the 
spectrum is continuous the possibility exists that a frequency lying midway be- 
tween two correlation frequencies would be missed. Therefore, a technique 
which would give a smooth approximation to the total input spectrum was sought. 
One such technique has been derived from adaptive filter theory. 

This project has been mainly concerned with the application of the adaptive 
filter technique. Briefly stated, this technique is as follows. Given a set of 
samples representing a constant amplitude signal plus noise, construct a sample 
data filter (composed of denominator terms only). This filter is to be such that 
if it were driven by white Gaussian noise, the output spectrum of the filter would 
give an approximation to the spectrum represented by the given samples. The 
sharpest pole of the filter should then represent the signal which was corrupted 
by the noise. Although the filter is derived mathematically by minimizing the 
mean square error of the output of the inverse filter when driven by the given 
samples, the implementation of the actual method is straightforward. 

The derivation of this method plus results from computer runs have been 
submitted to GSFC. 
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AN EXPLORATION ON HANSEN'S METHOD 
OF PARTIAL ANOMALIES 

by 

George E. McCluskey, Jr.^ 


Hansen's method of partial anomalies was investigated in order to de- 
termine if it would be suitable for use in the calculation of the perturbations of 
artificial satellites and space probes having highly eccentric orbits. This method 
has received very little attention since its publication in 1856. The main reason 
for this appears to be the lack of celestial bodies having highly eccentric orbits 
with the orbital parameters determined with sufficient accuracy to warrant a 
perturbation calculation. With the advent of the space age, this is no longer the 
case. 

The commonly used methods of celestial mechanics express the perturba- 
tions of the perturbed body in trigonometric series involving multiples of the 
mean anomaly, eccentric anomaly, or true anomaly. If the eccentricity of the 
perturbed body is less than about 0.67, the series will converge with sufficient 
rapidity, although for e > 0.5, the calculations become very tedious. If e > 0.67, 
the convergence is extremely slow or may even fail completely. Consequently, 
for large eccentricities purely numerical techniques must be used. 

Hansen’s method of partial anomalies introduces quantities known as partial 
anomalies. If the perturbations are expressed as trigonometric series involving 
the partial anomalies, the convergence may be made very rapid for any value of 
the eccentricity. The method also has the great advantage over the numerical 
methods of being semi-analytic. 

In consequence, Hansen's method of partial anomalies is well suited for the 
computation of perturbations of satellites, such as the AIMP and IMP satellites, 
and interplanetary probes with highly eccentric orbits. 

This work was published in X- 552-6 7-3 81. 
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MAXIMIZING DATA STORAGE IN A DIGITAL DELAY LINE 

by 

Dr. John B. Peatman^ 


For special purpose digital systems design a delay line can be used as a 
low cost memory, where the storage of 1000 - 20,000 bits are required. This 
upper limit of 20,000 bits represents the present state of the art and is de- 
termined by the dispersion of a signal between input and output of the delay line. 
This limitation is dictated by the step response of a delay line. In essence, the 
input is differentiated, delayed by T d seconds and differentiated again at the 
output. Thus a step input should come out as a doublet. Dispersion spreads this 
doublet out so that T seconds fall between its peaks. The maximum storage 
capacity is determined by the rate at which step inputs can be put into the delay 
line and still be distinguished at the output. This limits the storage capacity to 
T d /T bits. 

Manufacturers of delay lines expend great effort to reduce this dispersion 
so as to maximize storage capacity. In this way 20,000 bits of storage have been 
obtained. 

The purpose of this study is to determine the problems associated with in- 
creasing this storage capacity through signal processing. We desire to apply 
step inputs to the delay line at such a rate that the output waveforms overlap. 
Then we wish to put the output through a filter such that the effect of each input 
step is apparent at the output of the filter. Several possibilities exist: 

1. Inverse filter 

2. Matched filter 

3. Waveform compression filter 

The inverse filter would recreate a step at the filter output for each step 
input to the delay line. Consequently if F(s) e~ T ds is the transfer function of the 
delay line, then the transfer function of the inverse filter is G(s) - l/F (s) . 
There are several drawbacks to the use of this filter. First, it is physically 
unrealizable since it must approach infinite gain as frequency approaches infinity 
in order to compensate for the attenuation of high frequencies inherent in F (s) . 
In addition, the inverse filter will include two integrations, which will cause 
baseline drift in the output in any real implementation of the filter. Finally in 
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the process of working toward an unattainable step output in response to a step 
input, any mismatch between G(s) and l/F(s) will be accentuated in an uncon- 
trolled fashion. That is, we have not designed the filter so as to minimize the 
sensitivity to such a mismatch. 

The term "matched filter" comes from radar signal processing and denotes 
a filter whose impulse response is 


g(t) = Kf (D - t) 


where K and D are arbitrary constants, while f ( ) is the response of the delay 
line previously defined. This filter maximizes the height of the signal component 
output pulse with respect to the RMS noise output due to white noise in the cir- 
cuit. For this problem this approach does not seem appropriate since noise is 
not the fundamental problem. The delay line plus output circuitry is relatively 
noise-free and predictable as compared with a radar system. Consequently it 
seems that a matched filter is based upon optimizing the wrong criterion. 

A waveform compression filter generated f (2t) in response to f(t), as- 
suming it is set up to compress the step response by a factor of two. Such an 
approach seems as if it is putting emphasis where it really counts for this 
problem. However, if it is assumed that f (t) consists of one period of a sinu- 
soid (a usefully simple idealization) then g(t), the impulse response of the filter, 
consists of a periodic waveform with complex shape. Such a waveform sounds 
unrealizable without getting into sensitivity problems. However since it is not 
necessary that the output be exactly equal to f (2t), and since f (t) itself is not 
really a finite waveform, the required g(t) may really be a decaying waveform. 
This seems as if it would really be necessary if the sensitivity to mismatch 
between f (t) and g(t) is not to be infinite. 



OPTICAL PROPERTIES OF LIGHT BEAMS 


by 

Harvey N. Rexroad ^ 

OPTICAL RESONATORS (THE CONFINEMENT OF LIGHT BEAMS) 

My NASA colleague, Nelson McAvoy, and I had agreed prior to my coming 
to Goddard that a good way to begin my tenure as a ASEE-NASA Fellow would 
be to review the available literature on the subject of optical resonators. Be- 
sides affording a means of introducing the author to the word of laser technology, 
there were two additional motives for the project. The first was education. 
Several ideas from physics and mathematics at the advanced graduate level are 
required for an understanding of this topic, which is of fundamental importance 
for laser research. Consequently, there is a need for an organized collection of 
all of this material. The second motive was directly related to the work in 
progress in the laboratory at Goddard. It was hoped that this study would provide 
an answer to questions having practical significance (the optimum geometry for 
maximum degree of coherence, the field pattern radiated from a coupling hole, 
the interpretation of data from a confocal scanning interferometer, etc.). 

A series of lectures and discussions, moderated by the author, did serve 
for the mutual educational benefit of everyone in the group. Notes have been 
completed and distributed covering the following topics: 

I. Introduction 
n. Basic Considerations 

A. Uncertainty Principal and Confinement 

B. Quality Factor 

IH. Introduction to Methods of Computing Properties of Optical Resonators 

A. Beam Iteration; Phase Transformers 

B. Self Consistent Field Condition 

C. Resonator Mode Condition; Diffraction Loss 
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IV. The Self Consistent Field Equations of Optical Resonators 

A. The Basic Diffraction Integral: 

Fresnel-Huggen’s Diffraction Integral for Optical Resonators 

B. Derivation of the Basic Diffraction Integral Using the Integral Rep- 
resentation of 15D. 

C. Simplified Basic Diffraction Integral. 

D. Identical Mirrors; Single Transit Field Consistent Equations 

E. Complete Pass Self Consistent Equations 

F. Symmetry Properties and Equivalence Relations for Optical Reso- 
nators 

Appendix A. Review of Electromagnetic Theory 

Appendix B. Gauss’ and Green’s Theorems: Delta Function Properties of 

V 2 (1/r) and Q 2 (1/r 2 ) 


Appendix C. Determination of Electromagnetic Field From Sources 

Appendix D. Kirchoff’s Surface Integral Representation 

Appendix E. Introduction to Integral Equations; Example of Solution of a 
Homogeneous Fredholm Equation of the Second Kind. 

Appendix F. Spheroidal Wave Functions 

The following additional notes will complete the coverage of lecture material 
and will be added in the near future: 

V. Solution of the Confocal Square Mirror Problem 

A. The Exact Solution 

B. Gauss - Hermite Functions 
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C. Field Pattern of Various Modes 

D. Propagation of Gaussian Waves 

VI. The Waveguide Approach 

VII. Image and Propagation Laws 

A. Image Rules of Fresnel Diffraction 

B. Propagation Laws 

C. Ray Optics 

D. The A BCD Law of Gaussian Imaging 

E. Circle Diagrams; Slide Rule Solution of Gaussian Beam Propagation 
Problems 

VIII. Practical Design of Optical Resonators 

A. Stability Condition 

B. Frequency Condition 

C. Mirror Spot Size 

D. Filling Factor 


IX. Perturbations 

Appendix G. Integration of Diffraction Integral for Beam Propagation 
Appendix H. Solution of Vector Wave Equation in Cylindrical Coordinates. 
Appendix I. Fresnel Diffraction Image Rule 
Appendix J. The Ray Transfer Matrix 

In Progress: A computer program is being written to determine the radiation 
zone field patterns from a non-confocal optical resonator with a hole for output 
coupling aperture. 
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SPREADING OF A LIGHT BEAM (Possible Publication) 


A somewhat natural extension of some aspects of the main problem led to 
the investigation of the spreading of a light beam as it propagates in space, 
which is analogous to the spreading of a quantum mechanical wave packet in 
time. The exact problem considered utilized classical electrodynamics to de- 
termine the Fraunhofer (far-field) diffraction due to a wave launched from a 
rectangular aperture at 0 in Figure 1. A gain for the finite surface, S, at F in 
the far field was defined as: 

G - (time average power through S per unit solid angle) divided by (time 
average power through A over 4-7t). 



F 



Area = S 
RECEIVE 


Figure 1 

I believe I have proved that a field distribution: 


E *„. (x - y) aS °» (c ' x/a) s - <c ’ y/t> exp fif 

with n - m = 0, over the aperture at 0 produces the maximum possible gain, G. 

In this equation S on (c,t) is a prolate spheroidal angle function, and the solid 
angle subtended by S to 0 is 4 c 2 A. 2 /(rr 2 A) . This maximum gain turns out to be: 


[R (i) 

LR oo 


(c,l)l 


where R on ( 1 ) is a prolate radial function of the first kind. This result agrees 
with a previously known one that G 0 = 4 -ttA/X 2 is the maximum gain when c = o. 
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There is, of course, a direct application of this maximum gain theorem to 
space communications and satellite tracking. It is also a fortunate circumstance 
that the lowest mode of a laser automatically produces this desired field distri- 
bution. 

Present Stage: A rough draft of a paper has been written. I am checking the 
literature to see if I have missed anything. Also, since the difference between 
this and previously known results is not large, there is some problem in assess- 
ing the value of this discovery. 
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TOWARDS OPTICAL COMMUNICATIONS SYSTEMS 

by 

Dr. Frederick Rojak^ 


Work during this summer was divided in roughly four categories. 

(1) Continuation of a study of three-dimensional optical data processing 
using holography. The advantages of holography over conventional 
photography became apparent for the case where the third dimension 
is traded for the acquisition of a large amount of two-dimensional in- 
formation. 

Attention was given to the problem of recording the parameters of 
spatially modulated coherent wavefronts, 

(2) An explanation of the theory of variable scale optical correlators. 

(3) An attempt at contributing to the development of an optical parity gen- 
erator for a convolutional encoder. 

(4) Giving two lectures on coherence theory and working toward a mean- 
ingful definition of '’coherence length." Familiarized myself somewhat 
with Glanber's work on this. 
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REAL TIME CONTROL SYSTEMS 


by 

Linus Sc hr age t 


As part of the development program for an automated, computer controlled, 
satellite tracking station a study was made in order to outline the functional re- 
quirements of the various components of the system. A substantial effort was 
devoted to the specification of a language for programming a real-time control 
computer. The proposed language has a standard alogrithmic language such as 
ALGOL or FORTRAN as a proper subset but has additional statement types 
particularly useful in real-time control programming. Examples of these state- 
ment types are the CONNECT statement which allows one to specify that a given 
external interrupt will cause a specified sub-routine to be entered when the 
interrupt is triggered. There are also statements to assign priorities to sub- 
routines and statements to disarm and disable interrupt channels. The execu- 
tion of subroutine at a future point in time can be scheduled through the use of 
the CAUSE statement. A measure of multi-programming is possible through 
use of group sub-routine calls and the DELAY statement. Utility subroutines 
can be compiled into re-entrant or pure code by the use of the PURE PROCEDURE 
declaration. 
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